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Abstract — We present a complete formalism for final-state (timelike) dipole-antenna 
showers including fermion masses, but neglecting polarization and finite-width effects. 
We make several comparisons of tree-level expansions of this shower algorithm to 
fixed-order matrix elements for hadronic Z decays, up to and including Z — > 6 par- 
tons, to which the algorithm can be consistently matched over all of phase space. We 
also compare to analytical resummations at the NLL level. The shower algorithm has 
been implemented in the publicly available VINCIA plugin to the PYTHIA 8 event gen- 
erator, which enables us to compare to experimental data at the fully hadronized level. 
We therefore also include comparisons to selected observables in 6-tagged Z decays. 

1 Introduction 

The large phase space opened up by the LHC is rekindling interest in the collider phenomenology of 
heavy coloured particles. Appreciable samples of top quarks with large Lorentz boosts are becoming 
accessible for the first time; energetic top and bottom quarks are sought as decay products of new- 
physics or Higgs particles; and coloured new-physics particles may also themselves give off radiation, 
though at suppressed rates close to threshold. 

Indeed, for massive particles produced near threshold, most of the radiation produced results from 
the violent deceleration of the incoming (massless) colour charges, i.e initial state radiation is dominant. 
The effects of multiple soft emissions from the massive partons themselves are then largely unimportant 
for the description of the event as a whole and only become relevant to define precisely the mass of the 
produced particle [1-3]. 

However, for the production of boosted heavy quarks or other coloured particles, either directly, via 
decay, or through gluon splitting processes, multiple emissions cannot be neglected. Mass corrections 
will generate differences in the shape of the evolving jet, and in the energy loss of the evolving particle. 
These effects must be taken systematically into account if we are to rely on physics models of these 
phenomena to distinguish "signal" from "background" production sources. 

On the theory side, calculations of observables involving massive particles present a unique set of 
challenges. The introduction of an additional scale in the problem, for each non-zero mass, leads to an 
increased number of terms in amplitudes, to modifications to the pole structure caused by the massive 
propagators and to more complicated phase-space boundaries and kinematics. The presence of massive 
final state particles shrinks the size of the phase space available for additional QCD emissions, both in 
fixed-order calculations and in parton showers. 
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The modifications to the pole structure imply different infrared limits in the massive case. In particu- 
lar, QCD radiation from massive particles can lead to soft divergences but cannot lead to strict collinear 
divergences, since the mass is acting as an infrared regulator. Traditional Monte Carlo (MC) shower 
descriptions, which rely on the relative dominance of collinear-enhanced terms, therefore become in- 
trinsically less accurate when non-zero masses are involved. Though it is possible to systematically 
improve shower descriptions to take into account universal mass effects (as, e.g., in [4]), one would still 
expect a relatively larger uncertainty from non-universal and/or subleading terms than in the massless 
case, simply because the leading singular behaviour itself is less strong. Consequently, corrections from 
higher-order matrix elements, generally referred to as "matrix element matching", may be relatively 
more important. 

In this paper, we shall attempt to address a relevant subset of these challenges, in the specific con- 
text of matched time-like dipole-antenna showers [5,6]. We restrict our attention to unpolarized stable 
massive particles, deferring a detailed treatment of helicity dependence, as discussed recently by [7] 
and finite-width effects to a future study. Still, our approach has some advantages. To the best of 
our knowledge, this is the first time a rigorous and systematic approach to mass effects has been in- 
corporated in an antenna-based shower Monte Carlo code [5,8,9]. We also generalize the fixed-order 
antenna functions derived in [10-14] to include variations in their non-singular behaviour and extend 
the unitarity-based matching formalism presented in [ ] to include tree-level matrix elements with up 
to four additional massive partons beyond the Born level. Due to the unitary nature of the matching 
corrections, this prescription can be used also in the soft and (quasi-)collinear regions and hence we 
expect the subleading properties of the resulting shower to be improved. This is a feature which is not 
possible with other approaches to multileg matching, such as MENLOPS [ 5], CKKW [16, 17] and 
related approaches [18-20], or MLM (see [21] for a description). The speed of the resulting matched 
calculations is also greatly improved as compared to the existing approaches, as discussed in [22]. 

Corrections at the next-to-leading order level have not yet been included in this work. We there- 
fore do not attempt to distinguish rigorously between different possible mass definitions, such as "con- 
stituent" vs "pole" vs "running" masses [2,3,23]. For the purpose of our studies here, we treat parton 
masses simply as effective parameters, to be determined from data. It has been argued that this should 
be comparable to using a perturbative mass definition evaluated at a scale of the order of the infrared 
shower cutoff [ ], though the corresponding scheme is only defined numerically by the shower algo- 
rithm. We expect that NLO matching for massive fermions will be able to provide some further insight 
into this question, but that is beyond the scope of the work presented here. 

This paper is organized as follows. In section 2, we discuss the factorization, kinematics, and in- 
frared limits of a single 2—^3 splitting involving massive partons. In section 3, we introduce the 
additional ingredients required to turn this into a framework for parton showering, including a discus- 
sion of trial functions and veto algorithm steps. The generalization of our evolution variables to the 
massive case and the treatment of g — > qq splittings in the shower are also addressed. Sections 4, 5, and 
6 then present comparisons to fixed-order matrix elements, to analytical resummations, and to 6-tagged 
experimental data, respectively. We round off with conclusions and an outlook in section 7. 

We note that the work reported has been made publicly available as a plug-in to the PYTHIA 8 event 
generator [25], starting from VINCIA version 1.026 [26]. 
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2 Massive Phase-Space Factorization and Massive Dipole-Antennae 



The dipole-antenna formalism [12,27-29] is constructed from two basic ingredients: 

1) an exact momentum-conserving and Lorentz-invariant phase-space factorization based on 2 — > 3 
mappings between on-shell partons 

2) a set of antenna functions that, combined with the phase-space factorization, capture the leading 
singular behaviour of gauge field theory amplitudes [28]. 

We return to how these are implemented in the shower context in section 3. In this section, we focus 
on a single "elementary" 2 — > 3 branching, e.g., as it would appear during a single step in a shower 
algorithm, and/or in the context of an antenna-based NLO calculation. We here focus on the gener- 
alizations necessary in the massive case, with details on the massless treatment available in [5, 6, 12]. 
We begin by giving some conventions concerning the notations we use, in section 2.1, then turn to the 
phase-space factorization in sections 2.2-2.3 and finally discuss the structure of the antenna functions 
in sections 2.4-2.5. 



2.1 Notation and conventions 

Given momenta p a ,Pb of massive particles a and b with masses m a and nib, it is convenient to use the 
notation, 

Sab = 2p a -p b = (p a + p b ) 2 - ml -mf , (1) 

which we adopt throughout this paper. With this notation, the relation expressing the conservation of 
the total centre-of-mass (CM) energy in a massive 2—^3 branching, IK — > ijk, becomes 

m2 iK = (Pi + Pk) 2 = sir + mj + m 2 K = + s jk + s ik + m 2 + m 2 + m 2 k . (2) 

The momenta involved in this branching are either called parent or pre -branching momenta for /, K and 
daughter or post-branching momenta for i, j, k. Hence we may express, the dot product of two daughter 
momenta i, k as, Si k = m 2 K ~ s ij ~ s jk ~ m i — m2 ~ m \- 

We shall also work with scaled invariants, normalized to the CM energy of the dipole-antenna given 
by m 2 IK 

Sab 

Vab = — o— > (3) 

mj K 

and scaled mass values, 

,a = ^. (4) 

m IK 

We denote dipole-antenna functions by the symbol a (for antenna) and represent the partons partic- 
ipating in the branching process IK — > ijk by subscripts ajjjx- The normalization of a is such that 
|M n+ i| 2 ~ a|M n | 2 in the relevant soft/collinear limits, to be elaborated upon in the sections below. We 
also define a corresponding colour- and coupling-stripped dipole-antenna function aj/ IK . The relation 
between a and a is 1 

a j/IK = ^Cj/IKO-j/IK (5) 



'Note that compared to the equivalent relation between these two antenna functions presented in [ ], the normalization 
to the phase-space factor is not present here anymore. We prefer to keep the phase-space factor normalization outside the 
dipole-antenna functions which are defined from matrix elements squared only. 
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where Cj/jx denotes the appropriate colour factor for the branching, as follows: Cf for qq — > qgq, Ca 
for gg — > ggg, either Cf or Ca for qg — > qgg, and Tr for the splitting of a gluon into a quark-antiquark 
pair. In the conventions for the colour factors we use [ J.wehaveC^ = Nc,Cp = 2Cp = Nq — I/Nq, 
Tr = 2Tr = 1 which makes the difference between Cf and Ca explicitly colour-subleading. 

Note furthermore that the lower index of the dipole-antenna function ol^/ik fully specifies the par- 
tons involved in the branching process, and hence those functions are uniquely determined. In particular, 
the radiators IK are also uniquely identified, and our antenna functions are therefore related to so-called 
"sub-antenna" functions in the context of fixed-order subtraction [12]. 



2.2 Phase-space factorization 

To define the dipole-antenna phase space characterizing the massive 2 — > 3 branching process /, K — > 
i,j, k, we consider the exact factorization of the n + 1-particle phase space d<J> n+1 into a n-particle 
phase space d$> n and a dipole-antenna phase space given by d$ 3 / d<I> 2 

i \ ak i \ d *3 (j?hPj,Pk) 

d $n+l (Pi • • • ,Pi,Pj,Pk, ■ ■■Pn+i;q) = d$ n (pi... ,Pi,PK, ■ ■■Pn+V,q) Afh , — r- • (6) 

d$ 2 (Pi,PK) 

In this equation, d$ n corresponds to the phase space for n outgoing particles with momenta p\ , ..p n and 
masses m\..m n , with total four-momentum q^ . In d<I> n only the parent momenta pj, pk appear. The 
relation between pi, px and pi, pj, pk, typically called "momentum mapping" in fixed-order subtraction 
contexts and "recoil strategy" in parton-shower ones, will be discussed in section 2.3 below. 

The dipole antenna phase space d$ 3 m is proportional to the three-particle phase space. It involves 
only the pre- and post-branching momenta pi,pk and Pi,Pj , Pk respectively and is given by [5, 6, 13], 



where the Kallen function A is given by 

A(a, b, c) = a 2 + b 2 + c 2 - 2(ab + bc + ac) , (8) 

(f> parametrizes rotations around the P/-Pft--axis in the centre-of-mass frame. As long as we restrict 
ourselves to unpolarized processes which we do in this paper, all emission probabilities are independent 
of (j). The factor d(j) / (2tt) will therefore be suppressed in the following. 

The boundaries of the three-particle phase space with general masses follow from momentum con- 
servation and the on-shell conditions. They are given by 

2m i m j = s~- < Sij < s+- = (m IK - m k f - m| - m 2 (9) 



s tkM~ 2( SlJ+ m 2 + m 2 ) 



[Sij + 



2m 2 ) ([sfj + 2m k (m IK - m k )] 



4 ~ ( s ij ) V 4 ~ Si i V 4 + Am lKm k - Sij 



(10) 
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Equivalently, these boundaries characterizing the physical phase space for the daughter partons i,j,k 
are determined by requiring the positivity of the Gram determinant A3 defined as, 



A I 22 2 2 2 2 I/ .22 2 N \ si n 

A3 = — [SijSikSjk - s i:j m k - s ik m j - s jk ra i + 4m i m J -ra fc ) . (11) 



2.3 Phase-space mappings 

To specify the phase-space factorization in equation (6), a momentum-conserving mapping, or "recoil 
strategy" in the parton shower language, that relates the three on-shell daughter momenta, pi, pj, and p k 
to the two on-shell parent momenta pj and px is needed. In a dipole- antenna approach, the radiators / 
and K can be both emitter or recoiler and the radiation emitted between them is shared smoothly and 
symmetrically amongst them. The mapping between daughter and parent momenta presented below will 
reflect this fundamental property. Independent of whether the momenta involved are massive or not, the 
on-shell and momentum-conserving 2—^3 mapping is not unique except on the boundaries of the phase 
space. Instead, there is a one-parameter family of such mappings. 

Generalizing the analysis in [28] to the case of non- vanishing particle masses we start by relating 
the momenta of the daughter particles to those of the parent momenta as follows, 

pi = xpi + rpj + zp k (12) 
Pk = (1 - x)pi + (1 - r)pj + (1 - z)p k (13) 

with the on-shell conditions 



22222222 22 
Pi = m i , Pj = m j , Pk = m k > Pk = m K , Pi = m i ■ ( 14 ) 

We use these on-shell conditions to re-express the parameters x and z in terms of the single free param- 
eter r and obtain 



2(4A 3 + m^(4-( S - )2 )) L 



(4-(^) 2 + 4Ai i ) 



+R ( sfj + 2m IK m k - Sij) + 8r (A 3 - mj K Aij] 



(15) 



' S 2 (4 - (s- k ) 2 + 4A iJfc ) 



2(4A 3 + mf if (4-W 2 )) 



-R (s^ k + 2m IK mi - s jk ^j + 8r (A 3 - m 2 IK A jk ) 



(16) 
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where we have defined 



R 2 



'IK 



A 



16A 3 [m IK 
2 mj mx 



2 1 2 
m IK + mj 



(-1) • det 



(-1) • det 



r (l _ r ) - (1 - r)m| - rm^] + [4 - ( s " ) 2 ] [ S ] K - ( S J R ) 2 ] (17) 

(18) 
(19) 



ml 





S 3fc 


2 




&ik 

2 


m 2 


m 2 


2 




s jk 


2 


2 



1 

4 



2sijm 



(20) 



(21) 



These equations characterize our 1 -parameter family of massive mappings. The parameters z and 
x are related to each other with the replacements i «->• k and R — > — i?. Contrary to the massless case, 
however, i? 2 > 0, (which corresponds to real momentum fractions x and z), is not true for arbitrary 
values of the momentum fraction r. 

We note that the massive dipole mapping of [ ] corresponding to a dipole made of the partons i, j 
and k which play respectively the roles of emitter (i), emittee (J) and spectator (k), is obtained as a 
special case, by setting r = x in the above formula. In this case we have, 



2m 2 IK 



+ 



( 



s-j + 2m IK m k 



'IK 



'IK) 



2mjj^\ 



\jK m2 IKi(Pi+Pj) 2 ,™>l) 

where, in the rest frame of pi + pj, we can rewrite the Kallen function using 

K m ]Ki (Pi+Pj] 



(22) 



\2 2\ 



AiPi+ptfEfal 



(23) 



where E k denotes the energy of the parton k, Vfc its three-velocity. As in the massless case, the massive 
dipole mapping is asymmetric under the interchange of the particles i and k, but symmetric under 
the interchange of i and j. While it is appropriate to use this mapping for a shower based on CS- 
dipoles (CS:Catani-Seymour) which distinguishes between emitter i, emittee j and spectator k, it would 
clearly be inappropriate to use it in a dipole-antenna shower like VINCIA where the roles of % and k are 
interchangeable. The dipole mapping is mentioned for comparison only. 

To get a geometrical picture of the mapping used in VINCIA, it is convenient to express the free 
parameter r in the massive mapping family presented above in terms of the angle ip between the daughter 
parton pi and the parent parton pi (see, e.g., [5, 8]). To this end, we write down the 4-product pi • pj as: 



Pi - Pi 



EiEi 



\pi\\pi\ cos ip 

1 



4(4A 3 + m 2 x (4-( S -)2) 



£ 2 (4 - (^) 2 ) ( 



s: k + 2m IK mi 



+8r [ Sj k + 2miKmi - Sjk ) A 



(s; k ) 2 + 4A ife ) 



(24) 
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Looking at equations (17) and equations (24), we see that if we approach the boundaries of the phase 
space (for example if we consider a soft emission pj — > or if we take the (quasi-)collinear limit for pi 
and pj), the Gram-determinant A3 tends to zero and the dependence of pi • pi on the free parameter r 
drops out. 2 

Inside the parton shower VINCIA, to define appropriately the 2 — > 3 branching, we need to fix the 
mapping. In other words, we need to fix the functional form of the free parameter r. If all particles are 
massless, the default mapping used is given by [28] 

r = . (25) 

Sij ~r Sjk 

This mapping has the properties that the interchange i k corresponds to r o (1 — r) and the 
momentum fractions are restricted to x > 1, < r < 1 and z < 0, where r = corresponds to the 
collinear limit Pj\\pk and r = 1 corresponds to the collinear limit Pi\\pj. 
In the massive case, we shall consider the following mapping, 



+ 



S jk s jk 



S V S ij + S jk s jk 



S 2 r 



l y 'IK ~ ( S Ik) 2 Sjk ~ S jk - (Sjj - Sij ) ^ 

2mj K 2m| K Sy — + Sjk — s- k 



where the condition 



> <r<r + = - 2 (27) 

2mj K ^ m iK 

ensures that R 2 > 0. 

Apart from reducing to the massless mapping in equation (25) for vanishing masses, this phase-space 
mapping has the "swapping" property that i k combined with I K corresponds to r -R- (1 — r). 
For mi = rrii and rriK = mk it also satisfies x > 1 and z < and can therefore be viewed as a 
generalization of the massless mapping given in equation (25). 

In figure 1 , the mapping given in equation (26) is illustrated in a Dalitz plot of the three-particle 
phase space of the daughter momenta pi , pj , p^. The phase-space boundary is marked with a solid grey 
line. Insets show the orientation of the daughter momenta for a branching with Sij and Sjk given by the 
centre of the inset, in the CM frame of the parent partons, with is chosen such that the radiated particle 
is moving "upwards". A mass configuration characteristic for Qq — > Qgq, with uiq = 0.25 uiik, 
m q = (left), is compared to the massless case (right). Notice that the physically allowed phase space 
shrinks considerably in the massive case, and that the invariant sij can only vanish in the soft limit 
Pj — > 0. The limit — > with j hard is not accessible. The mass effects on the mapping are most 
pronounced for configurations which are close to the edge of the phase space and far away from the soft 
limit. For the rest of the phase space they are relatively unimportant. A similar illustration for massless 
partons can be found in [6] . 



' The sole exception to this occurs for Sit = s ik where the angle tp does depend on the functional form of r. 
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(a) Qq — > Qgq (Q massive, q massless) 



(b) qq — > qgq (q massless) 



Figure 1: Dalitz plot of the dipole-antenna phase space for IK — > ijk for massive partons (left) as compared 
to massless ones (right), using the scaled invariants y a ^ defined in equation (3) as coordinates. The boundary of 
the physically allowed phase space is drawn as a solid grey line. Insets show the orientation of the ijk momenta 
corresponding to the centre of each inset, in the CM frame of the parent partons, with parents oriented horizontally 
and 4> chosen such that the gluon is radiated "upwards". The mass values used in the left-hand pane are to/ = 
rrii — 0.25 wiik , ttik = rrik — rrij = 0. 

2.4 Pole structure 

Since masses act as infrared regulators in the collinear region, the pole structure of massive amplitudes 
is actually simpler (less divergent) than that of their massless counterparts. A specific example of this 
is given in figure 2, in which we show the ratio of the amplitudes squared for the processes Z — > QgQ 
relative to Z — > qgq, as a function of the Qg opening angle, for Mz = 91GeV, E g = lOGeV and 
niQ = 4.8 GeV (Q stands for a massive quark while q stands for a massless one). The dip in the thick 
solid line for 0y — > is generated by the mass-shielding of the collinear enhancements, relative to the 
massless case (thin line). 

However, the calculation of observables with massive final state particles still involves the treatment 
of potentially large mass-dependent logarithmic terms. They correspond to collinear divergences which 
are regulated by the quark mass, therefore they become divergent in the massless limit. For observables 
that are infrared safe in the massless limit, these logarithmic terms cancel in the final result, but they 
can still appear at intermediate steps of the calculation, for example in the separate evaluation of real 
and virtual contributions. They are of the form \n{Q 2 /m 2 ), where m is the parton mass and Q is a 
characteristic scale of the hard-scattering process. These mass-dependent logarithmic terms are related 
to the quasi-collinear [3 1 ] limit of the matrix element, the definition of which we shall recall below. 

In a fixed-order approach, the potentially large logarithmic contributions induced by mass terms are 
taken care of in the context of subtraction methods [32]; terms which mimic the singular behaviour 
of real matrix elements are added and subtracted. The construction of these terms relies heavily on 
the factorization properties of amplitudes in their soft and (quasi-)collinear limits [31]. In the antenna 
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Figure 2: Illustration of the dampening of the collinear singularity for Z — > QgQ: squared matrix elements with 
(thick) and without (thin) mass corrections, normalized to the massless case, as a function of the opening angle 
between the quark and the gluon, for constant E g = lOGeV and toq = 4.8 GeV. 



framework presented in [13, 14] (and in the dipole formalism [ ] that predates it), the main building 
blocks, massive antenna (dipole) functions and phase-space factorizations, are therefore constructed 
so as to reproduce exactly the quasi-collinear and soft behaviours of real radiation matrix-elements in 
the corresponding limits. For cross sections which are well-behaved in the massless limit, the explicit 
cancellations of the \n(Q 2 /m 2 ) -terms also ensure numerical stability in the limit m — > 0. 

For some observables which are not infrared safe in their massless limit, such as ones sensitive to 
the details of the fragmentation process for example, the cancellation of the mass-dependent logarithms 
is incomplete. Terms of the form tie ln n (Q 2 /m 2 ) appear in every order of the expansion. In the case 
of a large hierarchy m <^ Q, these terms jeopardize the convergence of the perturbative series. It is 
necessary to resum them to all orders to obtain a meaningful result, as is done, for example, for the b- 
quark fragmentation process in [33], to which we compare the massive VINCIA dipole-antenna shower 
in section 5. However, in order to construct this shower, we must first consider the soft and quasi- 
collinear limits more carefully and define how the massless splitting functions and soft Eikonal factors 
are generalized in the presence of massive particles. 

The infrared singularity properties of tree-level colour-ordered matrix elements involving only mass- 
less partons have been well studied in [ ]. In the limit where a gluon j is soft with respect to its 
neighbouring partons i and k, the colour-ordered matrix-elements squared |A^ n +i| 2 for (n + 1) partons 
factorizes into a universal soft Eikonal factor Sijk and a colour-ordered tree-level squared amplitude 
where gluon j has been removed. For the squared amplitudes we have, 

\M n+ i(l,--- ,i,j,k,--- ,n + l)| 2 g 2 s Cij k S ijk \M n (l, ■ ■ ■ ,i,k,--- ,n + l)| 2 (28) 

where g 2 = 4Tra s is the strong coupling, Cijk is a colour factor that tends to Nc in the leading-colour 
limit, and the massless Eikonal factor is given by 

Sijk = ■ (29) 

Similarly when two neighbouring gluons or a quark and a gluon become collinear the colour-ordered 
matrix elements factorize. Depending on the nature of the partons involved different collinear factors 
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are obtained. Partons which are not colour-connected do not lead to singular behaviours of the colour 
ordered matrix-elements squared, hence the soft or collinear factors only involve the neighbouring par- 
ticles to which the unresolved particle is colour-connected. 

In the massive case, essentially the same factorization properties still hold, provided the collinear 
limit is generalized to the quasi-collinear limit (see below). For the emission of a soft gluon from 
massive radiators, the factorization of the matrix element into a soft Eikonal factor times a reduced 
matrix element with the soft gluon omitted works in the same way as for massless partons. The soft 
Eikonal factor given in equation (29) needs however to be generalized. Written in terms of the parent 
partem masses mj and and the invariants between the daughter partons i, j and k, the massive soft 
Eikonal factor reads 

o / \ 2s ih 2mf 2m\ 

S ijk {mi,m K ) = 2 2~~ ( 3 °) 

SijSjk s- k 

which has two new mass-dependent terms compared to the massless Eikonal factor defined above. 

The quasi-collinear limit of a massive parton with momentum decaying into two massive partons 
j and k is given by, 

p£->zp",p£->(l (3D 



with the constraints, 
at fixed ratios, 



P 2 = m 2 m . (32) 
Pj • Pk,mj,m k ,m jk -> (33) 



m 2 ml tn 2 k 
__3_ ^Z*_. (34) 

Pj ■ Pk Pj • Pk Pj ■ Pk 

The key difference between the massless collinear limit and the quasi-collinear limit is given by 
the constraint that the on-shell masses squared have to be kept of the same order as the invariant mass 
(Pj + Pfc) 2 > w i tn the latter becoming small. In these corresponding quasi-collinear limits, the colour- 
ordered [m + l)-parton matrix element squared factorizes into a reduced m-parton matrix element 
squared multiplied by quasi-collinear splitting functions, the latter are generalizations of the Altarelli- 
Parisi splitting functions [34] from which they differ by mass-dependent terms. In four dimensions, they 
read 



p , ._l + (l-z) 2 2m 2 



'J 



z s 



<!<:) 



, 2 (35) 



Pqq->G(Z, m q , S qq ) = Z 2 + (1 - zf — ^—^ . 

Sqq + 2m\ 

We now turn to a description of the full massive dipole- antenna functions as implemented in VINCIA. 



2.5 Massive dipole-antenna functions 

In general, the full forms of the dipole-antenna functions are obtained by normalizing a three-parton 
tree-level matrix-element squared to a corresponding two-parton squared matrix element, stripped of all 
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couplings and colour factors and normalized to reproduce the known collinear splitting functions and 
soft Eikonal factors in the corresponding unresolved limits. 

In the fixed-order context, specific sets of such dipole-antenna functions have been derived for the 
massless case in [10-12] and for the massive one in [13, 14]. In principle, there is an infinite set of 
similar dipole-antenna functions, differing by non-singular ("finite") terms and hence having the same 
soft and (quasi-)collinear limits, which could equally well be used to construct the subtraction terms. 
For example, consider gluon emission off a qq antenna. To cover the limiting behaviour of this emission 
process, we could use the matrix element for 7* — > QgQ normalized to the one for 7* — > QQ, 

fj-r^QgQ ( m 2 m \ _ 2s ik 2m l 2m \ 1 ( Sjj s jk \ 

Alternatively, we could use the process H — > QgQ, which gives 

H-+QgQ f m 2__ „.. _ m _\ _ 2s ik 2m l 2m l , 1 ( Sjj , S jk | 

$ / i 



~ a 9 /m i m ^' ^' s ^ m ^ m «) ~ S^7 k ~ ~gf s\ + SlK -2m a m, 1 + 



s ij s jk stj s* k sik - zm q rriq \s jk 

(37) 

In both of these expressions, the denominator factor sik + xm q mq is proportional to the two-parton ma- 
trix element to which the three -particle matrix element is normalized. These two different constructions 
would give rise to slightly different integrated and unintegrated subtraction terms, but the final result 
would in either case be completely independent of which one is used. 

For a parton shower, however, the behaviour of the dipole-antenna functions away from the phase- 
space boundaries is important to determine the amount of radiation produced. The most obvious ex- 
ample is that of adding a positive constant to a dipole-antenna function. This would result in a slightly 
higher rate for hard emissions in the parton shower (and consequently smaller Sudakov factors) with- 
out changing the limiting behaviour of the dipole-antenna function. In the context of shower uncertainty 
evaluations, it is therefore useful to generalize the definition of the dipole-antennae, to allow for continu- 
ous variations of the ambiguous non-singular terms, as done for the massless case in [ ]. The possibility 
of varying finite parts in the parton-shower framework is a particular and important feature of the VIN- 
CIA code. Other parton showers, whose evolution equations are based on fixed kernels, such as the 
Altarelli-Parisi splitting functions [34] and/or the Catani-Seymour dipole ones [35], do not provide this 
particular uncertainty measure. 

However, the presence of quark masses greatly increases the number of possible finite terms that 
could be added, hence we shall still place some limitations on the type of terms we will allow for, as 
will be described in detail below. As a starting point, we require that we must be able to reproduce 
the dipole-antenna functions which were derived from physical matrix elements, such as those given 
above. We then choose a generalization of the resulting parametrization, in such a way that the finite 
parts of all the dipole-antenna functions are parametrized in a similar way. We consider it preferable to 
have a rather general parametrization of the finite parts of the antenna functions because a change in the 
parametrization itself would require a change in the program code, whereas changes to individual terms 
within a given parametrization can be made without even recompiling the code. 

Subsequently, we must decide which values to assign the finite coefficients by default. Since we 
chiefly intend to use them for variations, our philosophy is to set most of them to zero from the start, 
allowing only for a few non-zero values to bring the tree-level expansion of the resulting parton shower 
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into reasonable agreement with the fixed-order matrix elements for Z decay up to Z — > 6 partons. Note 
that explicit comparisons to such matrix elements are given in section 4. 

In the context of our shower model, one must also require that the dipole-antenna functions be 
positive definite, since they act as branching probability densities. This is the case for all dipole-antennae 
considered in this paper. 

With the default choices fixed, we also define two antenna-function variations which we consider 
reasonably extremal, which we call "MIN" and "MAX". Our approach here has been to choose the coef- 
ficients for the MIN set as small as possible without introducing negative values for the dipole-antennae 
and then choosing the MAX coefficients such that the difference between the default coefficients and 
the MIN coefficients is at least as big as that between the MAX coefficients and the default coefficients. 
Note that we have not varied all possible finite coefficients in the MIN and MAX sets, but only a small 
subset of them, so it is conceivable that some physical variations could fall outside the range we define 
here. As always, uncertainty evaluations are more of an art than an exact science. We expect to learn 
more about the reasonableness of our choices as we expand to more processes in the future and can 
make explicit comparisons to more matrix elements. 

In the following we shall present the decomposition of the dipole-antenna functions used in VINCIA 
into their singular and finite parts. The following antennae are needed: gluon emission from a qq, qg, 
and gg parent antenna, and gluon splitting from a qg or gg parent. As mentioned in section 2.1, the hard 
radiators are always uniquely determined, and hence the antennae we discuss here are equivalent to the 
ones referred to as "sub-antennae" in fixed-order contexts (see, e.g., [12]). For each of these five antenna 
types, we give four separate sets of "finite terms": DEF (default), MIN, MAX, and GGG, with the latter 
reproducing the fixed-order antennae defined in [12-14]. Condensed summaries of the corresponding 
finite-term values are given in tables 1, for gluon emission antennae, and in tables 2, for gluon-splitting 
ones. 

To define the dipole-antenna for gluon emission off a massive QQ pair (where the quark and the 
antiquark may or may not be of identical flavour), we start with the generic form 

0>g/qq{ m IK> s jki rn 'q^ m q) = 

1 f ^Vik _ 2/4 _ 2/4 1 / yjj ^ y jk 

m 2 IK \yijyjk yl y) k 1 - n 2 q - n\ + x iM q Hq \yjk Vij 

where F g i qq represents an arbitrary "finite" function, i.e. a function which is regular in all soft and 
(quasi-)collinear limits. The antenna function derived from Z decay, as in the fixed order context and 
called j4g = Og corresponds to F g j qq = with x = 4 and m q = m q as listed in table 1 . 
We allow for the following optional terms in F g / qq , 

Fg/qq {yij,Vjk, fJ-q, M?) = Coo + C w ( yi j + Ujh) + C 2 o{yf j + y] k ) + C n yijy jk 

+ n q )(M$ + MlSiVij + yjk)) 
+04 + ^_)(M 2 ° + M™( yij + y jk )) 

+fi qN {M^ + Mil (Vij + Vjk)) (39) 

with the default values C a b = = 0. Note that this form of F explicitly respects charge conjugation 
symmetry (i k). In principle, one could allow for terms with higher powers of masses and/or in- 
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qq -»■ qgq qg -> qgg 



a g/qq 


Def 
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MIN 


MAX 


a g/qq 




Def 
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Cio 
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-1 
-1 


7.5 
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-7.5 
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Table 1: Colour factors and finite parts for four different examples of the colour-ordered gluon-emission an- 
tenna functions. The VINCIA default antenna set (Def) is compared to the "GGG", MIN, and MAX variations. 
Coefficients which are not listed (or which are represented by "-") are zero. 



variants, but for the simple purpose of uncertainty estimates, we believe the form above gives sufficient 
flexibility. 3 

The default QQ antenna function used in VINCIA is thus, 

3 (def) (m 2 \ _ i ( 2/4 1 ( yji , yjk\\ 

V« (rn IK ,s v ,s 3k ,m q ,m q ) - ^ ^— yf , y% + 1 _ ^ _ ^ { yjk + y J) 

(40) 

which corresponds to choosing zero values for the finite part and for x in equation (38). In the left-hand 
pane of table 1, we compare the default values in VINCIA to the GGG functions and to a "MIN" and 
"MAX" variation that we use for uncertainty estimates. Note that a non-zero M coefficient is introduced 
in the "MIN" case, in order to avoid negative regions in the massive dipole-antenna function. 

The influence of quark masses on the default QQ antenna function is illustrated in figure 3. This 
figure shows contours of constant values for the antenna a g / q q in a Dalitz plot of the three-particle phase 
space for massless quarks (short dashed) and for massive quarks (dashed). For massive quarks, the 
contour lines start to avoid the boundaries of the phase space (drawn as a solid grey line) for high values 
of the antenna function. This is a direct consequence of the presence of terms of the form {—sfj/m 2 ) in 
the massive soft Eikonal factor given in equation (30). 



3 The denominator factor sik + x m q mq in equation (38) has x = in VINCIA, but we retain the possibility to vary it for 
uncertainty estimates. 
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3'ij 



QQ -> QgQ 



Figure 3: Dalitz plot showing contours of the massive (dashed) and massless (short dashed) gluon-emission 
dipole-antenna function a g / q q, with tuq = 0.15 itijk- Contours of constant values of the dipole-antennae are 
shown for a = 5, 50. Grey solid lines denote the boundaries of the massless and of the massive phase space 
respectively. 



For gluon emission in a Qg dipole-antenna, we use the generic 4 form 



l 9/l9 { m IK7 s iji s jki m q) 



1 



2j4 1 

4 



m IK \ DijVjk 

We allow for the following finite terms, 

V?5 



Vjk 



Vij 



+ 



2 



(41) 



= Cqo + C w yij + C 01 y jk + C 2 oyfj + C 02 y] k + CnVijl/jk 



+/i g (M 1 + Mi y„ + M \y jk ) 
+fi 2 q (M 2 + Afg, + Mq 2 ! + M 2 



10 + 

yij 

+ M3 i _ u M^ + m4 i _ 11 ^) . 



(42) 



Hi k Vij Uik UijUik 

In a fixed order context, the corresponding massive was derived from the matrix element for the decay 
of a neutralino into a massive gluino and two gluons, \ 999 [1 3], using an effective Lagrangian [10]. 



4 There is an additional ambiguity which originates from the collinear gluon singularity. In VINCIA, we use identified 
particles and therefore we have to distribute the collinear gluon splitting singularity pj || onto two dipole-antennae, one in 
which j is a hard radiator and only k can become soft and one in which k is a hard radiator and only j can become soft. When 
doing so, we could introduce aj f> fc-asymmetric term of the form x{y*j ~ yik)/{yjk(l ~ fJ%)] to the dipole-antenna function 
which contributes in the collinear limit j || k. To first order, such an asymmetric term cancels out, but it would still influence 
the shower at higher orders. However, since values different from \ — tend to lead to negative dipole-antenna functions for 
high quark masses, we have so far not enabled the option to vary it in the VINCIA code. 
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qg -> qq'q' 
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gg gqq 



a q/gg 


Def 


GGG 
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MAX 
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JU 00 
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Table 2: Colour factors and finite parts for four different examples of the colour-ordered qg —> qq'q' and gg —> 
gqq antenna functions. The VINCIA default antenna set (Def) is compared to the "GGG", "MIN", and "MAX" 
variations. Coefficients which are not listed (or which are represented by "-") are zero. 



Similarly to the massive dipole- antenna a g / q q, the denominator factor s 2 K = (m 2 K — m 2 ) 2 is pro- 
portional to the coupling-stripped two-particle matrix element for neutralino decay into a gluino and a 
gluon, x —> 99- It has mass dimension four in this case because the neutralino-gluino-gluon coupling 
has mass dimension —2. The parametrization in equation (42) contains finite terms which are propor- 
tional to Hg/d/ijUik)- Although these terms seem out of place in a quark-gluon dipole-antenna function, 
they are indeed part of the fixed-order antenna d® if the quark is massive. Their appearance is connected 
to the fact that the physical matrix element for \ 999 from which d® is extracted is symmetric under 
the interchange of the two gluons. We are using the quark-gluon dipole-antenna function to generate 
emissions in a situation which is decidedly asymmetric between the two gluons since one is identified 
as the hard radiator and the other is the emitted particle. For this reason, the terms proportionate to 
/j, 2 1 (yijyik) have been deactivated by default in VINCIA and are not even considered for the purpose of 
uncertainty estimates. 

The values for the other coefficients are summarized in the right-hand pane of table 1, with the 
coefficients reproducing the fixed-order d$ given in the "GGG" column. 

Obviously, there are no (leading-order) mass effects for gg — > ggg. We include the corresponding 
generic 5 form of the dipole-antenna function here for completeness, 

a g/gg (m 2 IK , Sij ,s ]k ) = 4- + M^lM + VJkO-V*) + ^ c a b \ (43) 

where the last term represents F g / gg , with coefficients as given in the bottom left-hand pane of table 1. 
The "GGG" column tabulates the coefficients of the f® function derived in [ ] using an effective 
Lagrangian for the process H — > ggg normalized to H — > gg. 

5 If a term xilJij ~ Vik) /[Vik{^ — fig)] is added to equation (41), we would also need to introduce a term xiVii ~yik)/yjk + 
XiVjk - Vik)/yij in equation (43) to ensure that a g / qg (pi,pj,p k ) + a g/aa (pj,p k ,pi) reproduces P g / gg forpj || p k . 



15 



We proceed with the dipole-antennae describing the gluon splitting processes. The gluon-splitting 
process Qg — > QQ'Q' is described by the following generic form, 



, 2 \ 1 1 f sf k + s 2 - 2m 2 , \ m 2 JK 
a q'/qg V m iK, s ih s jk ,m q , m ql ) = - — 2 + , + -j-Ff/qg , (44) 

z b jk T \ S /X J* q' J IK 

with finite terms, 

-^'A/s = 2 (^oo + Cio + Cbi yjfc + C 2 o y| + C 2 yf fe + C X1 yijyjfc) 

+ ; ^% ( M oo° + Mf ° y - + M \V + M? (M 2 ° + M 2 ° ^ + Afg?y ifc )) 



£t 2 

H — 9 2 (Coo + Cio + Cbi yjfc + C 20 y 2 ,- + C 02 yf fe + Cn y^y^ 



/! "' ''" ^Mq 1 " + Vij + + fi q (m 2 ° + M 2 ° Vij + M$y jk } ) 



„4 _ 

+ *' 2 (<? + M°o y^ + Mo^y^) , (45) 

Vjk + V ' 

whose values are listed in the left-hand pane of table 2. In a fixed order context, the massive e5j quark- 
gluon antenna can be derived from the decay of a neutralino into a gluino and quark-antiquark, i.e from 
the process \ —> GQQ- The resulting partons, gluino and a quark-antiquark pair can either be massless 
or massive. The presence of two different masses is the main reason why the structure of the finite term 
Fq' /qg i s more complicated than for the other dipole-antenna functions. The values of the coefficients 

4 

column of table 2. 



which reproduce the sub-antenna e% for the most general case (all partons massive) are given in the GGG 



Contours of constant value of this function with default parameters for the finite part are illustrated 
in figure 4, for three different mass combinations, shown from left to right. Dashed contours represent 
massive functions, short dashed contours represent massless functions. Apart from the modifications 
to the size of the physical phase space, the mass effects are only numerically important away from the 
collinear limit. 

Finally, the gluon-splitting process gg — > gQQ is described by the following generic form of its 
dipole-antenna function 



Ml 99 { m2 IKi s iji s jki m q) 



^ik 



+ s 



'J 



2s 



+ 2m\ 



+ 



2m 2 



IK 



s jk + 2m 2 



+ 



ni- 



ne 



*IK 



Fq 1 /qg 



(46) 
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(a) Qg Qq'q' 



(b) qg gQ'Q' 



(c) ->■ QQ'Q' 



Figure 4: Dalitz plot showing contours of the massive (dashed) and massless (short dashed) gluon-splitting an- 
tenna function a^i / qg , for three different combinations of massive and massless partons, with tuq = 0.15 tjiik- 
Contours are shown for a = 1, 4 in (a) and (c) and for a — 0.5, 2 in (b). Grey solid lines mark the boundary of the 
phase space. Contrary to the gluon emission case, there are no qualitative changes caused by the introduction of 
quark masses. 



with finite terms, 



— — 2 (Coo + Cio yij + C i yjk + C 20 yfj + C 2 y% + c n yijVjk) 



+ 



Mo 



Coo + Cio yij + Cqi yjk + C 2 o + C02 + C11 yjjyjfe 



+ 



%7c + 2^ 



(47) 



In this case, the corresponding fixed-order antenna can be derived from the decay of a Higgs into 
a gluon and a massive quark-antiquark pair i.e from H — > QQg. It can be obtained from the dipole- 
antenna analogue given in equation (46) by setting the finite parts F according to the right-hand pane of 
table 2. 

Note that the finite parts of the the two dipole-antennae related to gluon splitting a^, / qg and aq/ gg 
have been parametrized in the same way. In the aq/ gg case, some simplifications occur though due to 
the presence of a massless parton in the final state. 

Finally, let us mention that VINCIA is not necessarily restricted to describe processes in the Standard 
Model with massive fermions in the final state. It could equally be used to describe processes with 
massive final state particles with different spin-statistics properties. Since many models of physics 
beyond the Standard Model contain new heavy coloured particles which can be scalars, VINCIA needs 
a default dipole-antenna function for those. Since the soft Eikonal factor of equation (30) is spin- 
independent, it can be used as a default dipole-antenna in VINCIA for those cases: 



1 2 2 2 \ 

OEikonalV 771 /^) s ij> s jki m li m K) 



2s 



ik 



2m] 



$ij$jk 



2m 2 K 

s 2 

s jk 



(48) 
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where mi = rrii and juk = rrik are the masses of the radiating partons and Sik can be obtained from 
the other arguments using equation (2). The colour factor used is equal to Cp = %Cf if both parents 
are in the fundamental representation, Ca if both are adjoints, and options for anything in between for 
mixed-parent antennae, with the default being \{Cf + Ca)- 

3 The VINCIA Formalism with Massive Particles 

In this section we present the main ingredients of our dipole-antenna shower formalism implemented in 
the present version of VINCIA. This formalism was first derived in [5, 6] to describe QCD radiation off 
massless partons. We here generalize it to take quark mass effects into account. For completeness, some 
aspects which are carried over from the massless case without modification are also summarized. 6 

VINCIA is based on the dipole-antenna picture of QCD radiation [ ] : its fundamental evolution 
step is a Lorentz-invariant 2 — > 3 branching process IK — > ijk by which two on-shell parent partons 
(/ and K) are replaced by three on-shell daughter partons (i, j, k), conserving four-momentum exactly. 
Dipole-antenna functions and phase-space mappings were discussed in sections 2.3 and 2.5. 

To construct an explicit shower algorithm, one must furthermore introduce an evolution variable 
(a.k.a. shower ordering variable or resolution scale), Qe, suitably generalized to massive particles as 
will be discussed in section 3.1. Together with this variable one also needs to define a concrete itera- 
tive strategy for generating and accepting "trial branchings" according to the Sudakov form factors, as 
described in sections 3.2 - 3.4. 

Matching to fixed-order matrix elements is performed using a unitary matching scheme, which re- 
sults in unweighted events matched to full-colour tree-level matrix elements squared, as explained in [6]. 
We give a brief summary of the main points of this method in section 3.5. VINCIA also calculates un- 
certainty estimates for its predictions at a moderate speed penalty, those are summarized in section 3.6. 

3.1 Ordering 

In the dipole-antenna formalism, subsequent emissions are naturally ordered by the nesting of the on- 
shell 2 — > 3 phase spaces: as more and more emissions (or g — )• qq splittings) are added, each dipole- 
antenna will, on average, carry a progressively smaller fraction of the total original centre-of mass energy 
squared s. We refer to this as "no ordering", but it could equally well be called "phase-space-ordering", 
since the only constraint implied is that energy-momentum inside the nested 2—^3 phase spaces is 
conserved. 

As demonstrated in [6,37], however, and further elaborated on for the massive case in our section 4, a 
dipole-antenna shower without any additional constraints would produce far too much radiation outside 
the double-logarithmic limit, i.e in the so-called "hard region". This is essentially due to the fact that 
an "unordered" shower approximation can be viewed as a sum of independent dipoles, while full QCD, 
beyond the 2 — > 3 level, has a more complicated multipole structure, with in particular destructive 
interference produced by colour-coherence effects. A reasonable agreement with all-orders QCD can be 
restored by enforcing a strict ordering of the emissions in terms of some "evolution scale", Qe- This 
variable represents a measure of (inverse) formation time or characteristic wavelength. 

6 We encourage readers unfamiliar with shower formulations to consult [3, 36] for recent pedagogical reviews. 
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The VINCIA formalism can accommodate a whole class of evolution variables which differ in how 
they prioritize soft emissions relative to collinear ones. We refrain from characterizing all possible 
evolution variables here, referring instead to the original publication [ ] for details. 

The two most important evolution variables implemented in VINCIA are the following: 



(Type 1) Transverse Momentum 



Ql = Ap\ 4™ = 4 — J = - 4 

m IK m IK m IK 



SijS jk _ A (2 Pi -pj){2 P j ygfc) m^=o (m\ - m 2 )(m 2 jk - m 2 ) (49) 
(Type 2) Dipole Virtuality 



s~.o 9 def r-i - / \ m j — ^ • / 9 99 9\ 

Q E = m D = 2mm(Sij, Sjk) = 2rmn{mp — mf^m^i. — mi 



(50) 



*ij "H ' ""jk k> ' 

where we again emphasize that the notation is used for the dot product Sij = 2pi • pj, which differs 
from the invariant mass squared mfj = {p% +Pj) 2 when non-zero rest masses are involved. Note that the 
rightmost expressions in the above equations are appropriate only to gluon emission, for which rrij = 0. 

While the imposition of such ordering conditions can extend the agreement with higher-order QCD 
to a much larger region of phase space outside the double-logarithmic limit (see, e.g., [6, 37]), it does 
have a side effect, which is formally beyond LL: in general, there will be small corners of the n-particle 
phase space which are not accessible through any sequence of strongly ordered branchings [37, 38]. 
Those are called "dead zones" , which the shower does not populate at all. These zones correspond to 
regions of phase space that are classified by the ordering condition as having no LL contributions. It is 
therefore consistent to set these corresponding contributions to zero at the LL level. 

As shown in [6], it is possible to avoid dead zones without re-introducing the large overestimates 
present in the "no-ordering" scenario. This can be done by allowing unordered branchings to occur 
with a suppressed probability which does not affect the LL accuracy of the shower. This is technically 
achieved by starting from an unordered shower and, instead of applying the strong-ordering condition 
as a step function, apply a smooth damping factor instead. We label this as an "improved" ordering 
condition, Pi mp , 



©strong— ordering ^ -Pimp {Qe > Q e) — ^ 9 > (51) 



Qj 
Ql + Q 2 f 



where Qe is the evolution scale of the current n — > n + 1 branching. Qe is a measure of the scale 
of the n-particle configuration, defined as the minimal value of the evolution scale evaluated over all 
partons in the n-parton parent configuration. For branchings which are at a much lower scale than the 
last one, i.e. the strongly ordered limit Q\ -C Q 2 E , this factor is unity, whereas branchings with the 
opposite hierarchy Q\ 3> Q E are strongly suppressed. At the point Q\ = Q 2 E , the suppression factor 
in equation (5 1) is equal to 1/2. We refer to this as the "smooth ordering" condition. 

Note that, since the P liav factor is everywhere smaller than unity, it can be applied as a probabilistic 
veto, which we make use of in the technical implementation. 

3.2 The Evolution Algorithm 

Formally, we can define a unitary evolution operator S({p} n , Q st art, Qstop) which generates the pertur- 
bative radiation off an n-parton state {p} n between the two resolution scales Qstart and <5 st op in terms of 
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an iterative Markov chain 7 



•S{{p\ni Qstarti Qstop) — A({p} n , Qstartj Qstop) 




aj/IK^ [{P}n, Qstart, Q l - 



) < ^({?'}n+l >yfe ' QrestarU Qstop) , (52) 



where the first line represents the fraction of states that remain unchanged by the evolution (i.e., the 
exclusive ra-parton fraction at the resolution scale Q st0 p), while the second line includes all states that 
do evolve (i.e., the inclusive (n + l)-parton fraction). To maintain unitarity, the second line must be 
minus the derivative of the first. Analogously to other time-dependent decay problems, the solution to 
this differential equation is that the Sudakov factor, A, must be the exponential of minus the integrated 
tree-level branching probability, 



with phase-space measures and dipole-antenna functions as defined in section 2. 

The sum in equation (52) and the product in (53) run over all possible 2 — > 3 branchings IK — > ijk. 
The integrals in equations (52) and (54) are performed over the range Q st art > Q l ^ k > Qstop, with 
Qijk _ Q E (jpi,pj,p k ) the evolution variable. 

The starting scale, Q start, represents the "factorization scale" of the n-parton configuration and may 
be given either by the invariant mass of the evolving dipole-antenna, by the restart scale defined by a 
previous branching, or by some externally imposed scale, depending on the type of ordering criterion 
imposed on the shower evolution. 

In a standard shower application, Q s top represents the infrared shower cutoff, or hadronization scale, 
with a value of ~ 1 GeV. For simplicity, we have assumed here that the definition of Q sto p in terms of the 
post-branching momenta pi, pj and p^ is the same as that of Qemit ■ In practice, this is not necessarily 
the case, see the section on hadronization in [ ]. 

The presence of S^p} 1 ^^ 1 ^ , Q^taro Qstop) hi equation (52) generates the continued (iterated) 
evolution of the (n + l)-parton state after branching, with Qjestart normally taken to be equal to Q l ^ k , 
for a traditional so-called strongly-ordered shower. 

Now we turn to the algorithmic steps themselves. Since the 2—^3 branching phase space is three- 
dimensional (two independent Lorentz invariants, e.g. Sij and Sjk, and the azimuthal angle which 
determines the global orientation), three independent random numbers must be picked for each step of 
the algorithm. The first of these is the Q em it scale, distributed according 1 — A(Q &un , Q em it) with the 
"no-emission" probability A given in equation (53). Next, the second independent Lorentz invariant 
must be generated. It is generated according to the integrand in equation (54). From the two Lorentz 
invariants and the azimuthal angle <j) (which is chosen according to a flat probability density because the 

7 Strictly speaking, strongly ordered showers depend on the scale of the last branching, and are therefore not completely 
Markovian. In the context of VINCIA, the difference is only really relevant in the context of matching to matrix elements and 
can therefore be ignored for the discussion of the pure shower algorithm. 



A({p}n, Qstart, Qemit) = Yi eXP (~ ^3 /IK (Qstart, Qemit)) 

IK— >ijk 



(53) 



■Aj/lK(Qstart, Qemit) = / T^Tf? a j/IK (Pi,PjiPk) , 



(54) 
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integrand does not depend on it), the momenta pi, pj, and pk can then be constructed from the parent 
momenta pi and px , based on the chosen kinematics map and the relations given in section 2.3. 

The algebra and computational overhead involved in generating the three branching invariants can 
be simplified considerably by judicious use of the veto algorithm (see [3, 39, 40] for pedagogical re- 
views). First, so-called "trial branchings" are generated, using a simplified form of the integrand in 
equation (54). A percentage of these trials are then rejected, using 

p _ a j/iK {mjjcSij^j^jm}) 

-'accept — I 2 \ ' W-V 

Atrial \P^IK ' ' ^Jk ) 

to determine whether a given trial branching should be accepted or not. Symbolically {m} stands for 
all the masses of the parent and daughter partons, clj/jk is the desired integrand in equation (54), and 
atrial is the simplified trial function. The veto algorithm ensures that the final answer has no dependence 
(apart from the speed with which it is obtained) on the form of the trial function used, requiring only 
that this function is an overestimate of the correct integrand over all of phase space, so that the accept 
probability equation (55) does not exceed unity. 



3.3 Trial Gluon Emissions 

For gluon emission, the trial function used in VINCIA is based on the double-pole singular behaviour of 
the soft Eikonal factor, equation (30), which it coincides with in the soft limit and overestimates every- 
where else. Using the notation conventions adopted in section 2.1, it can be adapted straightforwardly 
from the massless case, and is given by 

_ &S r 2mj K 

"trial-emit — ~7~^A > 

where the overestimates as > as and C t riai-emit = Ca > Cj/IK can t> e use d to guarantee sufficient 
"headroom" for arbitrary coupling constants and colour factors, respectively. Since the mass correc- 
tions to the soft Eikonal factor and to the quasi-collinear splitting function are negative (as given in 
equations (30) and (35) and illustrated, e.g., in figure 2), this function is also guaranteed to be an over- 
estimate in the massive case. 

In figure 5, we attempt to give a more concrete impression of the suitability of this trial function, 
for abb — > bgb branching in a dipole- antenna of mass y/s = 91 GeV (left pane) and = 45.5 GeV (right 
pane), with m& = 4.8 GeV in both cases and a gluon energy of E g = 10 GeV. These values were selected 
so as to yield plots that can be compared directly to those in [4]. The x axes show the gluon emission 
angle in degrees, going from the collinear limit (zero angle) at origin to a 90-degree emission angle on 
the right-hand edge of the plots. The value of the gluon-emission trial function defined in equation (56) 
is shown as a thick solid black line. It is everywhere larger than all the other curves and hence represents 
an overestimate, as desired. The Eikonal factor defined in equation (30) is shown as a thin solid line 
(red). The default VINCIA qq antenna function defined in equation (40) is shown as a thick lighter 
(yellow) curve; it is slightly larger than the Eikonal factor, but is still everywhere smaller than the trial 
function. For completeness, antenna functions derived from two different LO matrix elements are also 
shown (for H° and Z° decay, shown with dots and dashes, respectively). Note that the matrix-element 
curves are closer to the default VINCIA antenna function than to the Eikonal factor, in particular in the 
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<J7= 91 GeV, m b = 4.8 GeV, E g = 10 GeV 



^7= 45.5 GeV, m b = 4.8 GeV, E g = 10 GeV 





Trial 




Eikonal 




VINCIA (unmatched) 




H->bgb 




Z^bgb 




"bg (degrees) 




Figure 5: Various levels of approximation for QQ — > QgQ, compared to the VINCIA trial and default (un- 
matched) shower functions. Left: the size of the dipole-antenna functions as a function of emission angle for a 
gluon energy of 10 GeV, for y/s = 91 GeV and mq = 4.8 GeV. Right: the same for y/s — 45.5 GeV. Lower panes 
show results normalized to the default VINCIA dipole-antenna function. 



zero-degree region. This is due to the fact that the default VINCIA antenna function not only reproduces 
the soft limit, but also the quasi-collinear limit of the full matrix elements. 

Note also that the trial function is closest to the physical antenna functions in the left-hand pane, 
where the values of the antenna functions are ten times larger than in the right-hand one (notice the factor 
10 difference in the y axis scales). The overall efficiency of the trial algorithm, which is dominated by 
the regions in which the trial function is large, is hence quite reasonable. 

The corresponding evolution integral for trial gluon emission, equation (54), was defined for mass- 
less partons in [6], in terms of the variables Qe and (, 



-4-trial-emit(Qstart) Qe 



C A 



y 2 

i start 



,2 

emit 







(57) 



We have suppressed a trivial integration over <fi and changed variables from Sy and Sjk to Q 2 E and an 
arbitrary (linearly independent) phase-space variable (, with \J\ the Jacobian of the transformation. For 
the two Qe definitions discussed in section 3.1, a convenient definition of £ is 

Sin 



c 



+ s jk 



(58) 
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For both of the evolution variables considered above, the inverse relations can be written 



an = CQ 2 , s jk = (1 - C)Q 2 , (59) 
with the definition of Q depending on the choice of evolution variable. For Qe = 2p±, 

Q 2 = Q e miK , (60) 

while for Qe = rno, 

Q 2 = . n ■ (61) 

mm(C, 1 - C) 

Using these relations, a set of (Qe,C) values can be translated unambiguously back to the original 
phase-space invariants (sij,Sjk)- We therefore emphasize that, although the £ variable plays a role 
analogous to the z fraction of traditional parton showers, it here serves merely to (re)parametrize phase 
space; there is no explicit dependence on this choice [6]. 

Since the massive phase space is contained within the massless one, the massless phase space gen- 
erator developed in [ ] can be recycled for massive momenta. In this case, the points which do not 
correspond to physical massive momenta are rejected. We do have to take into account, however, that 
the two-particle phase space d$2(pi, Pk)> which normalizes the the dipole-antenna branching phase 
space, depends on the masses of the parent partons I and K and is not anymore just given by the invari- 
ant mass of the parent partons as in the massless case. Instead, we have 



d$ 2 (pj,;pk) d$ 2 (pi,p K ) 



(62) 

mj=mj(=0 



where the Kallen function A is defined in equation (8) and the reduced masses are defined as \i = 
m/miK, with m 2 IK = (pi + Pk) 2 ■ In the massive case, the relevant evolution integral can therefore be 
written as the massless evolution integral modified by a massive phase space factor , 

Arial-enrit(wfjo Qstart, Qemit) | m/=mjf=0 
^•trial-emit — > ) = • (p3) 



We take this into account by applying the Kallen factor as a multiplicative pre-factor on the trial emission 
probabilities. With this replacement, the treatment of trial gluon emission derived in [ ] in the massless 
case, can be carried over to the massive one without further modifications. 



3.4 Trial Gluon Splittings 

For gluon- splitting antennae, i.e. those containing a g — > qq branching, we again make use of the 
leading singularity structure of the underlying process. This improves on the treatment in [6], in which 
equation (56) was used for both gluon emission and gluon splitting. Irrespective of whether a gluon splits 
into a massive or into a massless quark-antiquark pair, the process is characterised by an s-channel gluon 
propagator with a singularity structure given by the factor l/m^. The effect of non-zero quark masses 
is only to restrict the emission phase space in such a way that the m q q — > singularity present in the 
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massless case cannot be reached with on-shell quarks. This suggests that the "optimal" trial function to 
use for gluon splitting processes (massless and massive) is: 



atrial-split = -7—nfTji — g~ , (64) 
4-7T rriqq- 

where nj denotes the number of active flavours and we recall that Tr = 1 in our conventions. 

While different choices for the evolution variable for gluon emissions are implemented in VINCIA, 
the evolution variable for gluon splitting has been fixed to gluon virtuality, m 2 qS , as also advocated 
in [41]. This is based on the fact that, while gluon emission involves a sum over terms that have different 
soft and collinear limits, with different evolution variables assigning different "times" to each region, 
here there is only one singular structure, in m? qg , and hence we consider this choice to be relatively 
unambiguous. Further motivation is provided by a comparison between mass- and p± -ordering given in 
section 4.2. 

What remains ambiguous is then the details of how to "interleave" gluon emissions and gluon split- 
tings in the parton shower evolution. The VINCIA algorithm sketched in this section is based on gener- 
ating a branching scale for every possible 2 — > 3 branching and choosing the highest of those scales to 
determine which branching occurs next. We therefore need a way to determine whether a possible gluon 
splitting is at a higher scale than a possible gluon emission. Once we have defined how we compare a 
gluon splitting scale with a gluon emission scale, the interleaving of gluon emission and gluon splitting 
works exactly the same way as the interleaving of gluon emissions from different dipole-antennae. In 
order to compare the scales, we use the fact that all gluon emission evolution variables as well as the 
evolution variable for gluon splitting, m qq , range from for strictly soft/collinear branchings to mjK for 
the hardest branching process kinematically allowed. This suggests that we can compare m qq directly 
with the gluon emission variable, for example 2p± . There is an ambiguity however in how we define this 
comparison. One could for example equally well compare mjK{m qq /mjK) a (a > 0) with the gluon 
emission variable. In VINCIA, we have so far chosen to maintain the nominal evolution variable, with 
a = 1. 

Let us here consider the case where parton K is the splitting gluon (as in qg — > qq'q' or gg — > gqq), 
parton I can be either massless or massive. The direct equivalent of equation (57), the integrated trial 
antenna function for gluon splittings, is then 

•~ 1 /"Qstart \ 

Arial-split(W'/i<', Qstart, Qemit) = n f^R — 2 2 / ^ S 99 ^Sqg ft (65) 

171 IK ~ m I J Q e 2 mit 47r m qq 

where m 2 1K —m\ originates from the two-particle phase space volume (we have inserted tuk = m g = 0) 
and 1 / rrigg corresponds to the trial function for the gluon splitting branching. Since we fix Q 2 E = m 2 g , 
the most convenient definition for Q = £ sp iit is simply the other phase space invariant in the 3-particle 
phase space d$3, normalized by the mass of the mother dipole-antenna, 

o 

771 - 

Csplit = -f 1 , (66) 
m IK 

with the massless phase-space boundaries 

Cmin(yB) = , Cmax(y£;) = 1 ~ y/jjE ■ (67) 
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with, ue = Q%/™%k = m qq/ m 'iK f° r our s P ec i nc choice of evolution variable. We recall that the 
definition of £ has no physical significance in our formalism (neither for gluon emission nor for the 
gluon splitting process). 

As in the procedure for trial gluon emissions described in detail in [ ] and adapted for the massive 
case in section 3.3, we replace the upper limit on £ by an overestimate during trial generation, 

Cmax = Cmax(QiJ m in) | m =o ' 

such that the ( integral in a given "evolution window" (with lower boundary Qemin ) becomes simply a 
constant 

/"Cmax ^ [q2 

Wl) = dC = Cmax = 1 - (69) 

/0 V m IK 



The evolution integral, equation (65), then acquires the form 

1 . rQiaM dQ 2 E as 



Arial-splitC^/Jr, Qstart, Qemit) - nfT R - a Cmax / ~WT~~T~ (70) 



where /xj = ^f^- If «s does not depend on Qe, then the integrated trial function in equation (70) 
simplifies to 

Arial-split («lTl<r> Qstart) Qemit) = n fTR~ ft J^Cmax In [ ^f"" J (71) 

1 ~~ Mj 47r V V start / 

while if we consider a first-order running a s as a function of fi R = k^Q 2 E = ^m^-, then the integral 
becomes 

A ( n n \ 1 ; 1, ( 111 ( fc 'Qstart/A 2 ) \ 

^trial-spUt^.Qstart^emit) - ^ []^Q2^ ) ^ 

where the [ln(ln())] structure seen in equation (72) reflects the single logarithms generated by the 
antenna-function singularities folded with the logarithm coming from the running of a s . 

In our treatment of flavour thresholds, a heavy flavour is treated as active, i.e. it contributes to 
the running of as and is allowed to be created in gluon splittings as long as we have m qq - > tuq in 
gluon splittings and as long as Qe > rriQ - which is 2p± > mg in the default settings - in gluon 
emissions. The threshold m qq - = itiq chosen here instead of the kinematical threshold m qq - = 2mg is a 
consequence of our interleaving of gluon emissions and gluon splittings discussed above combined with 
the gluon emission threshold choice made in [ ]. Of course the kinematical conditions for the secondary 
production of the heavy flavour are always enforced. Therefore the only consequence of the fact that 
our flavour threshold is below the kinematical threshold is a slight loss of efficiency of the algorithm 
for rriQ < m qq - < 2m,Q due to vetoed gluon splittings into the kinematically disallowed heavy flavour 
region. The alternative of adjusting the gluon emission flavour thresholds such that the gluon splitting 
thresholds are at m qq = 2niQ - which would be equally valid at leading logarithmic accuracy - is not 
implemented at present in VINCIA. 

With this definition of A given in equation (70), the generation of trial branchings can be carried 
over from the formalism presented in [ ]. 
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3.5 Matrix-Element Corrections (Matching) 



The procedure for matching VINCIA to leading-order matrix elements [ ] is not affected by the presence 
of mass terms and can be adapted to the massive case by just upgrading the massless matrix elements and 
antenna functions to massive ones. In accordance with the antenna-factorization formalism, all particles 
are treated as being on shell, both in the antenna functions and in the matrix elements. 

Briefly summarized, the strategy is as follows. Similarly to the PYTHIA [ ] and GENEVA [43] 
approaches, the VINCIA matching formalism relies on the antenna shower itself to provide an all- 
orders phase-space generator that captures the leading behaviour of full QCD by construction. At each 
trial branching in the shower, the accept/reject probability can then be augmented by a multiplicative 
factor that goes to unity in the collinear and soft limit, but which modifies the branching probability 
outside those limits. The modification factor is constructed precisely such that the full-colour leading- 
order matrix element squared is obtained after summing over shower histories. The approach relies 
heavily on unitarity and is qualitatively different from other multi-leg approaches in the literature, such 
as the MLM (see [21] for a description) and CKKW [ ] ones. An important technical difference is that 
VINCIA only requires a Born-level phase-space generator, with all higher multiplicities being generated 
by the shower. There is therefore no need for separate phase-space generators for the higher-multiplicity 
matrix elements, which can result in significant speed gains, both in terms of initialization time (virtually 
zero in VINCIA), and in terms of running speed. The reader is referred to [ ] for further details. 



3.6 Uncertainty Estimates 

Another crucial point concerns how to estimate reliably the accuracy of the resulting calculation. Argu- 
ing that variations only of the renormalization scale is insufficient at best (and misleading at worst), a 
more comprehensive approach for all-orders (matched-shower) perturbative calculations was proposed 
in [ ] and implemented in VINCIA. As with the prescription for matrix-element matching, this approach 
can again be adapted to the massive case straightforwardly. 

Briefly summarized, VINCIA is able to compute a number of weights corresponding to alternative 
shower settings along with each event. The central weight, corresponding to the current user settings, is 
unity, while each of the alternative weights represents the relative probability that the event would have 
been produced when running with the corresponding alternative setting. The uncertainties are evaluated 
in a way that explicitly preserves unitarity, and hence the weights for a particular alternative setting 
average to 1 over a large number of events. 

The uncertainties accounted for in the present version of VINCIA do not differ from those presented 
for the original formulation [ ]. Those include, for each branching: variation of the renormalization 
scale by a factor of 2 in either direction, variation of the non-singular terms in the antenna functions 
from a "MIN" setting to a "MAX" setting (see section 2.5), variation of the shower evolution variable 
between p^-like and mass-like choices, and variations proportional to 1/Nq. In the context of matrix- 
element corrections (see above), variations of the numerical value of any "matching scale" applied can 
also be included. See [ , ] for further details. 
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4 Comparison to Fixed-Order Matrix Elements 



By construction, the massive dipole-antenna shower formalism presented in the preceding sections re- 
produces the (quasi-)collinear and soft limits of the amplitude squared for a single shower branching. 
In this section we present some examinations of its behaviour across multiple (combinations of) gluon 
emissions and/or gluon splittings. Specifically, we compare tree-level expansions of the shower to fixed- 
order matrix elements for Z — > 4, 5 and 6 partons, treating each (leading-)colour structure separately 8 . 

We consider three possible evolution orderings: no ordering, strong ordering in transverse momen- 
tum and smooth ordering in transverse momentum, as defined in section 3.1. A comparison with other 
orderings can be found in [ ] for massless partons. For gluon splitting, we also consider the difference 
between ordering in transverse momentum and ordering in gluon virtuality. The dipole-antenna func- 
tions are the default ones given in section 2.5, with the phase space mapping defined by equation (26). 
The matrix elements are obtained from MADGRAPH 4.4.26 [44]. 

For each parton multiplicity, we make a flat (uniform) scan over the relevant n-parton phase space 
using an implementation of the RAMBO algorithm [ ] provided in VINCIA. In each phase space point, 
the tree-level expansion of the shower weight, wps, is given by a sum over nested antenna functions, 
subjected to the selected ordering criterion. E.g., for Z — > qggq, the tree-level expansion of the shower 
weight is 



where Q e denotes the evolution variable and tilded variables are obtained by reclustering the final-state 
momenta to intermediate 3-parton states, using the inverse of the shower kinematics map described in 
section 2.3. The functions express the strong-ordering condition for each of the two possible clustering 
histories that lead from 2 to 4 partons in the shower. For an unordered shower, they would be absent 
(i.e., unity), whereas for a smoothly ordered shower, they would be replaced by the Pi mp factor defined 
in section 3.1. For higher numbers of partons, more terms are generated, for which we use an iterative 
code structure to compute the relevant sums. 

The tree-level expansion of the shower weight wps can then be compared to the norm squared of the 
appropriate fixed-order colour-ordered sub-amplitude squared, |.M^| 2 , forming the ratio 



This ratio thus represents an estimate of the relative accuracy of the shower (or, rather, its tree-level 
expansion) phase-space point by phase-space point. For simplicity, all couplings and colour factors are 
set to unity in this comparison. 

By studying how the distribution of R n evolves with n, we obtain a useful indication of the accuracy 
of the shower, and how this accuracy evolves with parton multiplicity. We consider these comparisons 
to be fairly conservative, since the soft- and collinear-enhanced regions only occupy a relatively small 
corner of phase space in a flat scan. Similar comparisons to tree-level matrix elements were carried 

8 Subleading-colour properties were studied for the massless case in [ ] and are not repeated here. 



^ps = [ag/ q g{q,gi,g2)a g / qq -{qgi,gig2,q)@{QE{qgi,gig2,q) 



Qe (9,51,52)) + 




(73) 
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out for massless showers in [6, 37, 38]. Here, we focus in particular on the modifications to these 
comparisons caused by non-vanishing masses. 

For comparison purposes, we will plot the logarithm of such ratios, log 10 (i?„), for different shower 
approximations, final state multiplicities, and parton masses. This logarithm gives a way of quantifying 
the amount of over or undercounting by the shower approximation. For phase space points for which the 
shower approximation reproduces the matrix element exactly, this logarithm is zero. For those points 
for which the logarithm is positive, the shower overestimates the matrix elements while it underestimate 
the matrix elements for negative values of this logarithm. 

4.1 Z -> QQ + gluons 

For massless quarks, the default qq — > qgq dipole-antenna function, which describes the first 2 — > 3 
branching in the shower evolution, coincides with the matrix element for Z — > qgq, and therefore the 
shower is "automatically" matched to the Z — > 3 matrix element. For massive quarks, however, the 
shower with default antennae differs slightly from the exact matrix element for Z decay already at this 
order. This was discussed in sections 2.5 and 3.3, with an illustration provided in figure 5. 

Turning to Z — > 4 partons and more, in figure 6, we illustrate how quark masses affect the dis- 
tribution of R n , for Z — > qq + 2 gluons (left), + 3 gluons (middle), and + 4 gluons (right), for three 
different evolution criteria: no ordering (top), strong ordering in p± (middle), and smooth ordering in 
P± (bottom). 

With no ordering (top row), the massless shower (solid black) has a large tail to the right, i.e., it 
substantially overcounts the matrix elements, as was also discussed in [6]. Towards the left of zero, it 
falls off extremely sharply (notice the logarithmic y axis), implying that the unordered shower approx- 
imation is an almost strict overestimate of the matrix elements. The introduction of masses changes 
this picture drastically, with a moderate ratio uiq/ E cm = 0.1 shown as a thin solid histogram and the 
larger rriQ/E cm = 0.3 shown with dashes: for both mass values, a lot of the overcounting on the right 
of zero is removed, and for the larger ratio the regions below zero, where the shower underestimates the 
matrix-elements, is populated. 

In the second row of figure 6, the introduction of strong ordering in p± improves systematically on 
the non-ordered approximation and the introduction of quark masses does not spoil this improvement. 
Even for quite large quark masses, the distributions remain almost centred around log 10 (i?) = 0. How- 
ever, as discussed in section 3.1, the price for strong ordering is the introduction of a dead zone, which 
we illustrate by plotting the underflow bin at log 10 (i?) = —2. Its size corresponds to a few percent of 
the phase space volume. 

In the last row, the change to a smooth ordering condition is illustrated. The dead zone is removed 
and, at least in the massless case, this smooth ordering condition further improves the agreement with 
the matrix elements relative to the strong-ordering case. It eliminates the tail of large undercounting 
that was present for Z — > 5 and 6 partons in the strong-ordering case and sharpens the peak around 
log 10 (-R) = also for Z — > 4. This improvement is much less significant in the massive case, but since 
the dead zone is still removed, we use the smooth ordering option as our default choice for massive 
partons as well. 

We note also that the effect of imposing strong ordering in p± is much more pronounced for massless 
quarks than for massive ones. To see this, we compare for instance the change in the black (massless) 
histogram between the top (unordered) and middle (strongly-ordered) left-hand panes of figure 6. We 
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Figure 6: Histograms of log 10 (R n ), as defined in the text, in a flat phase space scan, for n = 4 (left panes), n = 5 
(middle panes), and n = 6 (right panes), for not ordered (top row), strongly ordered (middle row), and smoothly 
ordered (bottom row) shower approximations. The parton shower uses the default dipole-antenna set defined in 
section 3. 
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investigate this further in the 2D phase-space plots presented in figures 7 and 8. 

Figure 7 shows the case for massless quarks. In the left-hand pane, no ordering condition is imposed; 
in the right-hand pane, strong ordering in p± . The axes of the figure have been chosen to be logarithmic 
in the two successive branching scales p± : i/E cm and p± t 2/p±,i, respectively, with p± i the emission 
scale of the first branching and p± 2 the emission scale of the second branching. Among the two shower 
histories for Z — > qggq, we show the p± values of the larger contribution on the plot. An average over 
the phase space points compatible with the corresponding values on the x- and y-axes is shown, obtained 
using the same flat scans of the phase space as for the one-dimensional phase space plots above. 

The effect of the strong-ordering condition is clearly visible in the right-hand pane of figure 7, 
removing the shower contributions in the upper half of the plot, corresponding to unordered branchings. 

Comparing the left-hand pane in figure 7 (massless case, no ordering) to that in figure 8 (massive 
case, with m/E cm = 0.3, no ordering), one notices that the part of the phase space in which the 
overestimate of the shower is biggest in the massless case corresponding to the region where the second 
emission is at a higher scale than the first one, is almost inaccessible in the massive case due to kinematic 
restrictions. 

In accordance with this, the effect of imposing strong ordering in p± is comparably small for heavy 
quarks, as shown in the right-hand pane of figure 8. We conclude that to impose an ordering condition 
is much more important for massless quarks than it is for massive ones. 



4.2 Including massless g — > q'q' splittings 

The gluon-splitting g — > qq dipole-antenna functions (a^i i qg and aq/ gg in the notation adopted here) 
only contain single poles in the region where the secondary quark-antiquark pair becomes collinear. 
Those antennae are therefore less singular than their gluon-emission counterparts, and hence there are 
intrinsically fewer g — ^ qq splittings than gluon emissions occurring (independently of the difference Ca 
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vs n/Tft in colour/flavour factors). However, for those gluon splittings that do occur, the smaller relative 
size of the (universal) singular terms, on which the shower approximations are based, as compared to 
possible non-singular (and non-universal) terms, imply that one can expect an overall worsening of 
the shower approximation for processes involving g — > qq, as compared to ones involving only gluon 
emission. An immediate consequence of this is, for instance, that the amount of strange and heavier 
quarks produced in perturbative g — > qq splittings is associated with substantial uncertainties in all 
current shower models. 

A first illustration of this feature is given by figure 9 which shows the ratio of the expanded parton 
shower (with no ordering condition imposed) to the tree-level matrix element for Z — > Qq'q'Q. The x 
axis now ranges from -4 to 4, rather than -2 to 2, allowing for a much larger range of shower-to-matrix- 
element ratios. This accommodates the most important feature in figure 9: the tail of high overestimates 
of the unordered shower approximation for a massless primary quark-antiquark pair has R values ex- 
tending up to approximately 10 4 as opposed to approximately 10 2 for gluon emission, for the same y 
range (i.e., same fraction of flat phase space). 

A phase-space scan similar to the ones shown in figures 7 and 8, revealed that most of the very high 
overestimates occur for configurations where the secondary quark-antiquark pair takes up almost all of 
the energy, which in turn forces the primary quarks to be soft. One possible shower history leading 
to such a particular configuration is obtained by having a collinear high-z gluon emission followed by 
a very hard g — > qq splitting, or represented pictorially illustrated in figure 10. Such occurrences are 
apparently all too frequent, in the unordered shower. Physically, we interpret this as a screening effect 
which is missing in the unordered approximation. By independently adding the splitting probabilities 
in each of the qg and gq antennae, we are not taking into account any screening effects produced by 
the collective qgq system, which become particularly relevant when two or more of those partons are 
collinear with respect to each other and hence should maximally screen each other. 

As illustrated by the histograms for light primary quarks (thin solid line) and for heavy primary 
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Figure 9: Histograms of log 10 (i?„), as defined in the text, for Z — > Qq'q'(gg)Q in a fiat phase space scan, for 
n = 4 (left pane), n = 5 (middle pane), and n — 6 (right pane). No ordering condition is imposed. 




quarks (dashed) in figure 9, the introduction of non-zero masses for the primary quarks improves this 
situation, since the configurations where the shower overestimates are the largest in the massless case 
simply cannot be reached for heavy primary quarks due to kinematical constraints. In terms of coher- 
ence, the strong dampening of the collinear singularity for massive emitters leads to an absence of the 
subsequent very strong coherence dampening that is present in the massless case. 

We now turn to ordered showers. By analogy with the case for gluon emission, figure 6, we expect 
that we can get rid of a significant part of the high shower overestimates for highly energetic secondary 
quark-antiquark pairs by imposing a strong ordering condition in the secondary quark-antiquark mass 
m q q for gluon splittings. An alternative choice that cannot a priori be excluded would be to use p±- 
ordering for gluon splitting as well. In figure 1 1, we make a first comparison of these two possibilities, 
for processes with primary massless qq parent partons and involving gluon splittings in a second branch- 
ing step. Either we use strong ordering in Qe = 2p± for gluon emission and ordering in Qe = m qq 
for gluon splitting (shown in the left-hand pane of the figure), or we use strong ordering in Qe = 2p± 
for both branching processes (shown in the right-hand pane). Specifically, in the left-hand pane, the 
subsequent gluon splitting is vetoed if m q i q i > 2p±(q,g, q), while on the right-hand pane, the gluon 
splitting is vetoed if p±(q,q',q') > p±{q,g,q) or p±(q,q',q') > P±{q,g,q) respectively. A sub- 
stantial over-counting for highly energetic secondary quark-antiquark pairs remains in the p± -ordered 
case (right-hand pane, top right corner), while an undercounting results when the evolution variable is 
changed to m qq in the second branching step (left). We conclude that the strong ordering in m qq does 
a better job of suppressing the high overestimates for the regions of the phase space where there is no 
leading-log contribution (top right of the plots). 

As in the case of gluon emission, we wish to avoid dead zones by switching to a smooth suppression 
of gluon splittings, using the suppression factor Pimp, defined in equation (51). However, in figure 12, 
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(a) Evolution variable m q q for gluon splitting. (b) Evolution variable p± for gluon splitting 



Figure 11: log 10 (wps/|A^| 2 ) for Z —> qq'q'q as a function of the secondary quark energies, rescaled to 
range from to 1. Left: Evolution variable m q q for gluon splitting, resulting in the strong ordering condition 
m q 'qi < 2p qq ^ qgq . Right: Evolution variable p± for gluon splittings, resulting in the strong ordering condition 
p q f~^ qq 9 < p qq ^ q9q . All quarks are massless. 



we illustrate that in this gluon splitting case, a naive application of this suppression factor -Pj mp , results 
in an overcounting for both choices of gluon-splitting variables. This figure shows the same distributions 
as in figure 11, but with the strong-ordering condition replaced by a smooth one. 

In a dipole-antenna shower that employs ordering in p± for all branchings, or for our smooth- 
ordering shower variant, an additional suppression mechanism is therefore needed to remove this over- 
counting and get reasonable agreement between approximated shower vs matrix elements for processes 
involving gluon splittings. The Lund dipole cascade implemented in the ARIADNE program [ ] uses 
the following factor to modify its gluon splitting antenna functions, 

_ 2m% 



mj K + m 



N 



where m/f K is the invariant mass squared of the parent antenna-dipole and m 2 N is that of the neighbour- 



ing dipole-antenna. Thus, if the preceding branching was collinear, with m 2 N — > 0, this factor produces 



a very strong suppression, while if the two dipole-antennae that share the splitting gluon have exactly 
equal sizes it goes to unity. 9 In figure 13, we show that the use of the PAri factor suppresses the overes- 
timates visible in figure 12 to a large degree, with a slightly better agreement obtained in the left-hand 
pane (for interleaved p± and m q q evolution) than in the right-hand one (with all processes ordered in 
PL)- 

The expansion of the resulting weights for Z — > 4, 5, and 6 partons with smooth ordering in 2p± for 
gluon emissions and in m q g for gluon splittings, are shown in figure 14a and for smooth ordering in p± 



'if the neighbouring dipole-antenna is much larger than the parent antenna, it even produces a slight enhancement, by up 
to a factor of 2, but this is more than compensated for by the reduction in the splitting probability of the neighbour itself. 
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(a) Evolution variable m q q for gluon splitting. 
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(b) Evolution variable p± for gluon splitting 



Figure 12: log 10 (wps/|A / i| 2 ) for Z — > qq'q'q as a function of the secondary quark energies, rescaled to range 
from to 1. All quarks are massless, smooth ordering is imposed, p± is used as the evolution variable for gluon 
emission. Left: Use m g g for gluon splitting. Right: Use p± for gluon splitting. 
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(a) Evolution variable m q q for gluon splitting. 

Figure 13: log 10 (i?4) for Z — > qq'q'q as a function of 
1. Smooth suppression of unordered emissions, with Pj 
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(b) Evolution variable p± for gluon splitting 
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; factor. 



for all branching processes in figure 14b. 

Although still far from the very good agreement obtained in the pure gluon-emission case (as com- 
pared with figure 6c), both the centre and the width of the weight distributions shown in figure 14b where 
smooth ordering including this additional suppression factor Pah is imposed for secondary massless qq 
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(a) Using p± for gluon emissions and m q q for secondary qq production. 
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(b) Using p± for both gluon emissions and for secondary qq production. 

Figure 14: Histograms of log 10 (i?„), for Z decay into a massive primary quark pair, a massless secondary quark 
pair and up to two gluons in a flat phase space scan, for n — 4 (left panes), n — 5 (middle panes), and n = 6 
(right panes). The P^i factor is used for gluon splittings. Top row: smooth ordering using p± as the scale for 
gluon emission and m q g for gluon splittings. Bottom row: smooth ordering using p± for both gluon emissions 
and gluon splittings. 



production, are now in tolerable agreement with the leading order (LO) matrix elements over a substan- 
tial fraction of phase space. Matching to the LO matrix elements can obviously be used to improve this 
agreement further, up to the orders for which matrix elements are available (see section 3.5). 

It should be emphasized, however, that the centre position of the distribution wps/ \-M\ 2 is highly 
sensitive to the finite parts of aqi/ qg . Since these pieces of ag// qg are not universal, the fact that the 
most frequent ratio of wps/ |7W | 2 is almost unity as demonstrated in figure 14a cannot be expected to 
be universal for all processes involving gluon splitting either. This is an unavoidable consequence of the 
less pronounced singular behaviour for these antennae, as compared to the gluon-emission ones. 

4.3 Including massive g — > Q'Q' splittings 

As a final set of comparisons we include massive g — > QQ splittings in the shower expansion. Since 
the corresponding dipole-antenna functions do not contain any poles at all (though one does appear for 
itiq — > 0), we should expect the shower approximation to be at its worst for this case, translating to very 
large uncertainties on, e.g., the amount of g — > bb splittings produced by it. 
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(a) Using p± for gluon emissions and m q q for gluon splittings. 
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(b) Using p± for gluon emissions and gluon splittings. 



Figure 15: Histograms of log 10 (i?„), as defined in the text, for Z decay into a massless primary quark-antiquark 
pair, a massive secondary quark-antiquark pair and up to two gluons in a flat phase space scan, for n = 4 (left 
panes), n = 5 (middle panes), and n — 6 (right panes). Unordered emissions are suppressed smoothly, for gluon 
splittings i-Ari is used. In the top row, the suppression factor for unordered branchings is calculated using p± as 
the scale for gluon emission and m q q as the scale for gluon splittings. In the bottom row, the suppression factor 
is calculated using p± for both gluon emissions and gluon splittings. Note: the parton shower uses the default 
antenna set. 



However, as illustrated by the plots in figure 15, the agreement is in fact at the same level as that 
obtained for massive parents in the previous subsection, except for a slight tilt of the distribution for 
rather heavy secondary quarks. We note that this is especially true for the "interleaved evolution" choice 
of using Qe = 2pj_ for gluon emissions and Qe = m-qq for gluon splittings, cf. figure 15 a, as compared 
to using p± for all branchings as illustrated in figure 15b. 

This strengthens our motivation for using the interleaved p±- and mass-ordered evolution as the 
default in VINCIA. Note also that we have checked that this conclusion appears to be robust against at 
least moderate variations of the antenna function finite terms. We conclude that there is still significant 
uncertainties surrounding massive g — > QQ splittings, but that the default choices made in VINCIA can 
at least be considered a sensible starting point. Of course, matching to matrix elements can still improve 
the situation, in particular for secondary quark-antiquark pairs of high invariant mass, by increasing the 
multiplicity at which the arbitrary finite parts of the antenna functions start to matter. 
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5 Comparison to Analytic Resummation 



Observables involving massive particles, like heavy-quark fragmentation processes, can be considered 
as being collinear-safe since collinear divergencies are regulated by the finite value of the heavy quark 
mass m. Thus, such processes can be computed order by order in perturbation theory. Nevertheless, as 
mentioned in section 2.4, mass-dependent logarithms of the form ln(Q 2 /m 2 ), where Q is the typical 
scale of the hard scattering process, appear at each order in perturbation theory. When the hard scale 
Q is much larger than m, these quasi-collinear logarithmic contributions can be large and have to be 
resummed to all orders to obtain reliable predictions for these observables. 

This resummation can be performed analytically, in a process-independent way, by using the per- 
turbative fragmentation formalism [46,47], which is summarized briefly below. Alternatively it can be 
performed numerically, using a parton-shower, such as VINCIA. 

In this section, our aim is to verify on a particular example, that the VINCIA implementation devel- 
oped in this paper, based on the exponentiation of massive dipole- antenna functions which reproduce 
the soft and the quasi-collinear limit of tree-level matrix-elements, performs the resummation of quasi- 
collinear logarithms correctly. We do this by comparing the predictions of VINCIA with those obtained 
from the analytic calculation of [ ] for the inclusive production of a single heavy meson H in e + e~ 
collisions. 

The single inclusive production of a heavy meson H as obtained by the process e + e~ — > V — > 
H{p)+X can be described by the production of a heavy quark pair Q — Q which subsequently hadronize 
to yield a heavy meson H. This hadronization process can be described by non-perturbative (heavy 
quark)-to-hadron fragmentation function (such as for example in the model of Peterson [48]) whose 
free parameters have to be extracted from the data. 

Since we are mainly interested in the perturbative contributions to the inclusive cross section, we 
shall here ignore this non-perturbative contribution. Consequently, we compare the predictions obtained 
from the calculation of the inclusive production of a heavy quark pair with our predictions obtained with 
VINCIA without hadronization. 

We consider the inclusive production of a heavy quark pair in e + e~ collisions in the kinematical 
region where the centre-of mass energy E cm of the collisions is much larger than the heavy quark mass 
m, i.e. E cm 3> m. At the same time the heavy quarks are produced in the perturbative regime with a 
mass m which is large enough so that (m 3> Aqcd) ■ We consider the following distribution, 

V(x,E 2 cm ,m 2 )^—^- (76) 

with o" tot the total hadronic cross section in e + e~. The energy fraction x of the heavy quark system with 
momentum p is given by x = 2 p ■ (p e + +p e -)/ E 2 m . 

Since we are mainly interested in the large x behaviour of the cross section, we concentrate on its 
so-called flavour-non-singlet contribution. In this case, the ingredients to the cross section are dependent 
on the difference of two flavour-non-identical parton species. 

In the following we shall first recall the main ingredients of the analytic calculation of [33] before 
presenting our comparison. 
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5.1 The heavy quark fragmentation formalism 

The inclusive production of heavy quarks in e + e~ annihilation is a purely perturbative process which can 
be described by the perturbative fragmentation formalism in which collinearly-enhanced contributions 
are resummed. In the large x region, logarithmically enhanced contributions due to soft radiation can 
also occur. A reliable theoretical prediction for this heavy quark production process can therefore only 
be obtained if both class of logarithms are appropriately resummed. This combined resummation has 
been performed in [ ] up to the next-to-leading logarithmic (NLL) accuracy. In the following we 
shall first describe how the collinear-enhanced logarithms are resummed before including soft-gluon 
resummation effects as well. 

Thanks to the factorization theorem of mass singularities the cross section for the production of a 
hadron H in e + e~ collisions can be written as a convolution of a process-dependent coefficient function 
C^ e and a parton-to-hadron fragmentation function denoted by D a / H , for each parton a = (q, q, g) 
involved, treating all flavours as massless. Performing a power series expansion in a s , one can compute 
each of the coefficient functions C^ e as a massless QCD partonic cross section and the fj,p evolution 
of D a /u, within perturbation theory. The fj,p evolution of D a / H is performed using the well-known 
Altarelli-Parisi evolution equations [34,49] which enable the resummation of logarithmic contributions 
of collinear origin in D a / H . Furthermore, provided D a / H is known at some initial factorization scale 
called hof,h, its form at any other higher factorization scale \xf,h can t> e determined from this evo- 
lution equation. Its form at the initial scale ^qf,h, chosen such that [j>of,h ~ Aqcd, is a purely non- 
perturbative contribution which has either to be modelled phenomenologically and/or determined from 
data. 

Within the framework of the perturbative fragmentation formalism, the cross section oq for the 
inclusive production in e + e~ collisions, of a heavy quark pair Q — Q is given as a generalization of 
the fragmentation formalism for hadrons described above, as follows: For each parton a = q,q, g, 
(still treating all flavours as massless), and for a given type of heavy quark Q (or Q), the heavy quark 
production cross section gq is given as a convolution of the previously defined coefficient function 
Cf^ e with D a jq, the perturbative fragmentation function of the massless parton a into Q. In the most 
general case, all massless partonic contributions (including gluonic ones) to coefficient functions and 
fragmentation functions have to be taken into account. As we restrict ourselves in this section to the 
non-singlet part of the cross section, only the primary production of a quark-antiquark pair of the given 
type will contribute, secondary production will be neglected. 

In calculations involving all-order resummation, it is usual to consider the equivalent expression in 
Mellin moment space where one uses the customary definition of the Mellin transform 

f N = f l &xx N - l f{x) (77) 
J o 

which transforms the convolution in x-space into a product in Mellin space. Using the perturbative 
fragmentation formalism as described above, the N moments of the inclusive distribution V may be 
written as, 

V {E 2 cm , m 2 ) = C^ +e ~ ] {a s (^),E 2 m , fi 2 , D N m 2 ) (78) 

where <7 to t/ ^ L °^ is given at order a s by (1 + a s /iT). e ^ denotes the Mellin transform of the e + e~ 
coefficient function while D n stands for the Mellin moments of the flavour non-singlet component of 



38 



the perturbative fragmentation function D(x), as denned in equation (77). Both of these contributions 
depend on the factorization scale fj, p which was introduced by factorizing the cross section. 

As mentioned before, we concentrate on the large-x behaviour, i.e x — > 1 (in momentum space) or 
equivalently, in Mellin space on the large- N limit (N S> 1) of this cross section. Note that, in [33] the 
inversion of the results to x space is performed using the Minimal Prescription of [50, 51], but we shall 
not use those results in what follows. We shall in fact compare the results obtained with VINCIA to the 
predictions of the analytic calculation obtained only in Mellin space, using equation (78). 

The resummation of collinear-enhanced logarithms of E 2 m /m 2 is achieved in equation (78) by writ- 
ing the N moments of the perturbative fragmentation function, denoted by Dn in this equation, as the 
product of an evolution operator En and D™, the perturbative initial condition for the heavy quark 
fragmentation function. Both of these functions depend on the factorization scale fip and on the starting 
point of the evolution called flop. The perturbative fragmentation function in Mellin space denoted by 
-Djv reads, 

D N {n 2 p,m 2 ) = E N (fjFp,/j,l F )D^(as(l4),l4,^op,m 2 ). (79) 

The evolution operator En is the solution of the Altarelli-Parisi evolution equations written in Mellin 
space as, 

din \j?p 

with the boundary condition En(\x^ f , l^p) = 1. In this equation, the anomalous dimension 7 is related 
to the well-known time-like Altarelli-Parisi splitting functions in Mellin space Pn [52-57], (up to the 
second order in a s ) by 

77V (a s) = gpr + (g) 2 pW + 0(4). (81) 

For heavy quark fragmentation, the starting point of the perturbative evolution fiQp is chosen to be of 
the same order as m, the mass of the heavy quark. As such it is a perturbative scale (i.e. m 3> Aqcd)- 
As a consequence, unlike the parton-to-hadron fragmentation function denned at the initial scale /j-of^h 
close to Aqcd , the initial condition D l ™ for the parton-to-(heavy quark) fragmentation function which 
depends on this perturbative starting scale fiQp, can be computed in perturbation theory as a power series 
in a s . At leading order, we have D mi (x) = 5(1 — x) (and therefore L>™ = 1), which expresses nothing 
but the fact that at leading order, a ^-flavoured jet is a b quark. As such the LO initial condition is trivial 
and independent of the mass m of the heavy-quark. 

In equation (78), the process dependence of the inclusive cross section is entirely contained in the 
coefficient function e . It denotes the Mellin transform of the e + e~ coefficient function in the 
MS scheme . It is given in [58]. For conciseness it is not presented here, its behaviour for large Mellin 
moment N is given by, 



■2 -In {oisinpj) E N (n F ,n 0F ) (80) 



#->-i + =^(-(^ +w -5) ta (f) + |^ + (! 



The perturbative fragmentation function Dn, as well as its components the evolution operator E^ 
and the initial fragmentation function D m \ are instead universal. Beyond leading order, those depend 
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on the factorization scales fip and /xof as well as on the renormalization scales /i and hq. We start by 
describing the required evolution operator En needed for our comparison. 

The resummation of large collinear logarithms is performed by solving the AP equations pertur- 
batively and by setting factorization and renormalization scales equal. This is realized by setting 
(ipsa fan E cm and /i « Hqf ~ m. 

Using in addition the second-order expansion of the anomalous dimensions as given in equation (81) 
the evolution operator for the non-singlet channel reads, 



p (0) 

En(» f , Mof) - {-^^ ) exp ^ (P N - —P N j j . (83) 

Using this equation, one sees that the leading collinear logarithms of the form (a s ln(E 2 m /m 2 )) n of 
equation (78) are resummed by combining the LO expression of the evolution operator with the LO 
expressions of the initial condition 7Jj n i and the LO expression of the coefficient function C e+e . Since 
at this order, these three ingredients do not depend on the heavy quark mass m, it corresponds to using 
a massless calculation to describe the heavy quark fragmentation process, which is a very crude approx- 
imation. The resummation of the next-to-leading collinear terms of the form a s (a s \n(E 2 m /m 2 )) n in 
the evolution operator given in equation (83) above requires the equivalent NLO expressions for these 
three components which are now mass-dependent. As a consequence, although the collinear-enhanced 
logarithms are formally resummed in VINCIA only at the LL level, to compare our predictions, we 
will use the full evolution operator En defined above in equation (83) which includes subleading loga- 
rithmic effects. The comparison can in fact only be made with the analytic result which shows a clear 
mass-dependent behaviour. 

To come to D N \ as mentioned before, the initial scale fiQp is perturbative as it is chosen of the order 
of the heavy-quark mass m and by this choice, the appearance of large logarithms of the ratio \iQpjm 
are avoided in the perturbative initial condition D N \ The perturbative initial condition up to order as 
which only resums the collinear-enhanced logarithms of the form E 2 m /m 2 is given below as, 

D^iasinl), fi 2 ,^ F ,m 2 ) = 



1+ 2, 



- C ? [ ^ " '> (£t ( ln (dffc? - l ) ) + ° W) ' <84) 



The corresponding expression in the large- N limit is given by, 



D N l {a s {no),li ,n Q F,m ) = 

1 + «S^o) Cp ^_ ln 2 N + _ + ^ lnN + Q (1) ^J + Q ^ (85) 

As can be seen from equation (82) for CV and equation (85) for D N l given in the large N limit, for 
1 — x <C 1 (corresponding to N S> 1), both the perturbative initial condition and the e + e~ coefficient 
function cannot be computed in fixed-order perturbation theory as they contain logarithmic contributions 
proportional to In N and In 2 N respectively. Those logarithmic contributions, arising through the radi- 
ation of soft gluons from the heavy quarks at the renormalization scale /xq ~ m, spoil the convergence 



40 



of the fixed order perturbative expansion at large N (or equivalently at large x) and have therefore to be 
resummed. 



Note that we expect the soft-gluon effect to be quantitatively more important for the initial condition 
of the heavy-quark fragmentation function Z?i n i than for the coefficient function C e+e since in the first 
case this soft-gluon or so-called Sudakov effect is controlled by the coupling o. s {^q) which is larger than 
a s (n 2 F ) present in the partonic cross section C e+e instead. 

In [33], this soft-gluon resummation was performed to NLL accuracy for both the coefficient func- 
tion and the perturbative initial condition. The expression of the coefficient function including soft-gluon 
resummation effects given in Mellin space, which we shall call soft-gluon-resummed part of the coeffi- 
cient function is denoted by Cfj reads, 



CUasifi 2 )^^,^,^) = exp (lnNgW(\) + g™ (\, |f , ^J)) (86) 



with A = 6o«s(/^ 2 ) hiiV and 



9 (1, (A) = ^(A + (l-A)ln(l-A)) 



» ,2) (*■ f ■ f ) - w (2A + 21,(1 " A) + ln2(1 " A)) + ^ ta<1 ~ A) (87 » 



1 M< 2 > EL\ A« E, 



(A + ln(l-A)) — -A^ln-^ -—Am 



2 



The coefficients and are given in the MS scheme by 

'67 7T 2 5 



aW = c f ,aW = \c f (c a ( 



18 6 9 1 



(88) 



The initial condition Dj n j including soft-gluon resummation which we shall call soft-gluon-resummed 
part of the initial condition is denoted by D™' S and reads, 



DZ*(a s (ti),ti,V Z oF,rn Z ) = exp ( InNgQ (\ ) + g£) ( A , ^, ^- ) ) (89) 

^0 MoF / - 



2 2 

m m 
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From the formulae above, one can see that the soft-gluon resummed parts of the coefficient func- 
tion and of the initial condition are not valid for arbitrarily large N values. Those have branch cuts in the 
complex N -plane, the former starting atN = exp(l/ (boas([i 2 ))), the latter atN = exp(l/ (26005(^0))) 
This behaviour signals the onset of non-perturbative physics at values of x very close to 1. Keeping this 
fact in mind, we will restrict our comparisons of the VINCIA results with the soft-gluon resummed 
calculation to moderate (N < 25) Mellin moments. 

The final step is to match the soft-gluon resummed part of the coefficient function, equations (86), 
and initial condition, equations (89), to the fixed-order result in such a way that the truncation of the 
matched result reproduces the fixed-order result and the logarithmic accuracy of the resummed part is 
preserved. This requirement does not fix the matched result uniquely, however. In [ ], the matching 
was performed additively. 

The final NLO +NLL resummed expression for the N moments of the (non-singlet) initial condition 
denoted by D™' matched is given by, 
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(91) 



In this equation, is given above in equation (5.1), -D™' S is the Sudakov-resummed part given in 
equation (89) while D 1 ™' S \ as denotes the expansion of this expression D™' S given in equation (89) 
to the first order in as- The definition of the matched e + e~ coefficient function is analogous to equa- 
tion (91), i.e. we have, 

^matched = ^ + _ (92) 

In this equation, Cn is the full e + e _ coefficient function at O (a s ) given in [ ], Cfj is the resummed 
part as given in equation (86) above and Cfj\ as corresponds to its expansion at order a s . 



The numerical difference of this matching scheme to the so-called "log-R matching scheme" used for example in [59] for 
the same observable is small for the parameter values we consider. 
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5.2 Comparison with VINCIA 



We now compare the analytical calculation of the fragmentation function V(E 2 m ,m 2 ), equation (78), 
described in the previous subsection, to VINCIA results transformed to Mellin space. For the analytical 
calculation, we use the evolution operator defined by equation (83) and consider two possibilities for the 
choice of coefficient function Cn and fragmentation-function initial condition D™. The simplest is to 
take the analytic results obtained using the fixed-order coefficient function given in [58], and the form 
for D™ 1 which only resums the collinear-enhanced logarithms given in equation (78). Alternatively, 
we use the equivalent matched expressions defined above in equations (91) and equation (92), which 
include soft-gluon resummation effects up to the NLL accuracy. 

In order to do the comparison appropriately, we need to fix a kinematical range where both the 
analytical results and the prediction from VINCIA are reliable. The range of validity of the analytic 
resummation calculation was discussed in the previous subsection. We expect it to work for a centre-of- 
mass energy E cm such that E cm 3> m and for m ^> Aqcd- For VINCIA, the requirements for the mass 
of the heavy quark m are given by, 

Qstop < m Q < E cm , (93) 

where <5 st op denotes the emission scale at which VINCIA would make the transition to a hadronization 
model. The first hierarchy Q s top "C mg is necessary to ensure that gluons which are soft compared to the 
quark mass can be emitted, whereas the second hierarchy rriQ -C E cm ensures that the quasi-collinear 
logarithms play an important role. 

Furthermore, since the accuracy of VINCIA is formally LL, there still remain considerable uncer- 
tainties. The uncertainty due to the choice of finite terms, estimated by using the "MIN" and "MAX" 
antenna sets defined in section 2.5, was found to be small in comparison to the uncertainty coming from 
the choice of renormalization scale. Hence we only consider the default antennae set for our results 
here. The renormalization-scale choice was then varied between = 2p± and fj,R = p±/2, to define 
the uncertainty range. (Note that, in the parton shower context, we do not have an explicit factorization 
scale, the evolution variable Qe has the role of a factorization scale instead.) 

The predictions obtained with the analytic resummation calculation will be taken at factorization 
and renormalization scale respectively given by fj,p = \i = Q and /io = fi-oF = Tn. 

The x-space fragmentation function obtained with VINCIA for a centre-of-mass energy of 10 4 GeV, 
quark mass vtiq/ E cm = 0.02 and as(m 2 z ) = 0.139 (using one-loop runing) is shown in the left-hand 
pane of figure 16, in the range 0.8 < x < 1. The right-hand pane shows the comparison to the analytical 
results in Mellin space, in the range 1 < N < 25, using the same parameters as for the VINCIA result. 
We find reasonable agreement between the two predictions over this range. The VINCIA prediction 
reproduces the main feature of the analytic resummation calculation. Due to its large uncertainty, it 
cannot however distinguish clearly between the different analytic predictions." We find comparable 
agreement for light (toq < 0.04 i^m) quarks, under the condition that the dependence of the VINCIA 
result on the cutoff C? s top is small. 



1 1 There are parameter values for which the resummation of soft-gluon logarithms has a much more pronounced effect than 
in figure 16. It turned out, however, that the dependence of the VINCIA results on the cutoff Q sto p is too large to allow for a 
comparison in these cases. 
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MR = nV±_ , 2p± 
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Band : VINCIA 1.026 (/j R = ±p ± , 2p ± , no hadr.) 

Solid Line : Analytic pQCD (incl NLL soft) 
Dotted Line : Analytic pQCD (without soft) 



Figure 16: Left: the T>(xb) distribution obtained with VINCIA with hadronization switched off, i.e., at Q stop — 
2p_Lstop = lGeV. The shaded band shows the uncertainty obtained by varying the renormalization scale in as 
by a factor of 2 and 1/2 respectively (with scale-stabilization switched off, see [ ]). Right: the Mellin transform, 
T>(N), including comparisons to analytic resummation with (solid) and without (dotted) soft-gluon resummation 
at NLL. E cm = 10 4 GeV, m Q = 0.02E cm , a s (m%) = 0.139. The VINCIA predictions use the default dipole- 
antenna set, strong ordering of the emissions in p±, no matching, and no secondary quark-antiquark production. 



6 Comparison to b-tagged experimental data 

As a final cross-check, we include three basic comparisons to data published by the SLD [ ], DELPHI 
[61], and L3 [62] experiments. We do not intend this to represent a full-fledged phenomenological study 
of b fragmentation in VINCIA. Rather, we wish to demonstrate that our implementation of mass effects 
in VINCIA yields sensible numbers also when compared directly to experimental data, including the 
effects of hadronization. 

We also include comparisons to default PYTHIA 8.150, on the same distributions. (This is trivial 
for us to do — we merely switch VINCIA off.) The resulting comparison represents a further validation 
and cross check (of both models) since the shower formalisms are quite different between PYTHIA and 
VINCIA, especially for massive particles. We also obtain a concrete check that the events generated by 
VINCIA are being treated by PYTHIA's string model of hadronization in a consistent manner. 

For PYTHIA, we use the default parameters of PYTHIA 8.150 [63]. We also note that default 
PYTHIA includes matching through Z — > 3 partons. For VINCIA, we use the VINCIA-specific tune of 
the light-flavour parameters reported in [ ] and include matching through Z — > 5 partons. Since tuning 
is not our main purpose here, we have not attempted to retune PYTHIA's 6-specific non-perturbative 
parameters to (re)optimize them for use with VINCIA. It is therefore possible that some further im- 
provements could be made in the non-perturbative description. 

In figure 17, we show the fragmentation function (a) for b quarks at the parton level 12 and (b) at 
the hadron level, compared to hadron-level SLD data for weakly decaying B mesons [60]. The top 
panes show the distributions normalized to unity, and the bottom panes show the ratios of theory to data. 
In the bottom panes, the inner (lighter) yellow bands indicate the statistical uncertainty, and the outer 

12 Obtained by switching off PYTHIA's hadronization model, i.e., the b quarks are evolved perturbatively to a scale Qstop = 
2p±had = lGeV (or 2pxevoi = 0.8 GeV in standalone PYTHIA). We ignore the slight difference in the exact value and 
definition of p± between PYTHIA and VINCIA, see [64] for a discussion. 
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Figure 17: VINCIA (thin lines) and PYTHIA 8 (thick lines), before (left) and after (right) hadronization, com- 
pared to SLD [60] data on the fragmentation function for weakly decaying B mesons 



(darker) bands indicate the combined statistical plus systematical uncertainty, added linearly. Note that 
the vertical error bars on the Monte Carlo predictions correspond to ±1.645a statistical uncertainty from 
the number of generated MC points, equivalent to 90% confidence. This is the default for how statistical 
MC uncertainties are displayed by VINCIA's plotting tool. 

Comparing the two panes of figure 17, we conclude that the non-perturbative corrections are signif- 
icant in the region above ig ~ 0.5, while the spectrum at lower xb is dominated by the perturbative 
prediction, for which the two codes are in good agreement, both with each other and with the data. In 
the high-x b region, VINCIA generates a slightly harder parton-level spectrum (peaked at higher x b) 
than PYTHIA. That is, VINCIA generates slightly less perturbative radiation. After hadronization, the 
VINCIA spectrum is somewhat softer than the PYTHIA one. The total non-perturbative component of 
the "parton-to-hadron" correction would therefore be evaluated as being slightly larger for VINCIA than 
for PYTHIA. The main properties of this correction are driven by the tuning to light flavours. As men- 
tioned above, a dedicated 6-specific tuning has not yet been carried out. Nonetheless, we note that the 
fragmentation spectrum obtained with the default tuning already appears to be in reasonable agreement 
with the data, cf. the right-hand pane of figure 17. 

In figure 18, we have used FASTJET [ ] to compare VINCIA and PYTHIA to DELPHI data 
[61] on , the ratio of the y23 distribution in 6-tagged events to light-flavour events, with 2/23 the 
dimensionless resolution scale at which the event goes from being a 2- to a 3-jet event, according to the 
kr clustering algorithm (using the E scheme). This distribution is sensitive to the value of the 6-quark 
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Figure 18: VINCIA (thin lines) and PYTHIA 8 (thick lines), before (left) and after (right) hadronization, com- 
pared to DELPHI [61] data on the ratio of j/23 in b- vs. light-flavour events. 



mass. Comparing left- (without hadronization) to right-hand (with hadronization) panes, we see that 
hadronization effects are large in the region below, roughly, 1/23 ~ 0.02. Above that value, the deviation 
from unity is generated mainly by perturbative quark mass effects. Both models use the default PYTHIA 
value of m& = 4.8 GeV, which appears to give a reasonable agreement with the measurement. 

In figure 19, we show a similar comparison for the Thrust event shape variable in 6-tagged events, 
compared to a measurement performed by the L3 collaboration [62]. Comparing (a) parton- to (b) 
hadron-level results, an interesting pattern can be seen, for both PYTHIA and VINCIA. The non- 
perturbative corrections are significant not only at low values ofr = 1 — T < 0.1, but also around 
r = 1/3, shown with a vertical dotted line in the upper panes of the plots. At this point, the distribution 
changes slope 13 ; thus, even though hadronization effects are parametrically strongly suppressed in that 
region, the small "smearing" they provide of the underlying perturbative prediction becomes relatively 
more important at exactly that point. With hadronization effects included, both models describe the data 
acceptably well. Combined with the good agreement found with the light-flavour variables included 
in [6], we conclude that the VINCIA code can at this point be considered validated as a Monte Carlo 
model for final-state showering and hadronization effects. 



7 Conclusions and Outlook 

A precise description of processes involving heavy coloured particles is of prime importance to the 
physics programme at current and future high-energy collider facilities. In this article, we have extended 
the timelike dipole-antenna shower formalism of [5, 6] to include massive fermions. Advantages of this 
treatment include the exact on-shell phase-space factorization which is inherent to the antenna formalism 

l3 Due to the underlying change from a 3-parton to a 4-parton quantity that the Thrust variable undergoes at that point. 
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Figure 19: VINCIA (thin lines) and PYTHIA 8 (thick lines), before (left) and after (right) hadronization, com- 
pared to L3 [62] data on Thrust in 6-tagged events. 



[12-14], the ability to vary the remaining ambiguous parts of the calculation in a similar manner as was 
already done for the massless case but now including explicit mass-dependent terms, and the smooth 
merging with the "GKS" matching formalism [6]. 

In this paper, we have shown evidence, through extensive comparisons to other calculational meth- 
ods, that the algorithm is physically sensible and that it can be expected to yield reasonably precise 
results over large parts of the phase space. Furthermore, using the GKS formalism, it can be matched to 
leading-order matrix elements over all of phase space up to any given fixed order (in practice, the limit 
is currently Born + 4 partons), which should help to give a systematic improvement even for fairly soft 
and/or collinear emissions. 

Although we have not attempted a dedicated tuning of the non-perturbative fragmentation parame- 
ters for heavy quarks, the preliminary tuning of the massless parameters reported in [ ], combined with 
the PYTHIA 8 defaults for the 6-specific ones, appears to give a reasonably good description of the B 
fragmentation function and of 6-tagged event shapes and jet rates at the Z pole. 

For a full-fledged application to physics at hadron colliders, two further ingredients remain to be 
developed: an extension of the formalism to spacelike (initial-state) showers, and the inclusion of finite- 
width decays, which can be important in the description of (chains of) resonance decays. Nonetheless, 
the step taken here is a necessary prerequisite, and can already be used to study questions involving 
final-state c and b quark fragmentation. 
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